Covariate-dependent Graphical Model Estimation via Neural Networks with Statistical Guarantees
Authors: Jiahe Lin
Yikai Zhang
George Michailidis
Venue: Accepted by TMLR ()

Abstract: Graphical models are widely used in diverse application domains to model the conditional dependencies amongst a collection of random variables. In this paper, we consider settings where the graph structure is covariate-dependent, and investigate a deep neural network-based approach to estimate it. The method allows for flexible functional dependency on the covariate, and fits the data reasonably well in the absence of a Gaussianity assumption. Theoretical results with PAC guarantees are established for the method, under assumptions commonly used in an Empirical Risk Minimization framework. The performance of the proposed method is evaluated on several synthetic data settings and benchmarked against existing approaches. The method is further illustrated on real datasets involving data from neuroscience and finance, respectively, and produces interpretable results.

Published in Transactions on Machine Learning Research (05/2025)
Covariate-dependenl Graphical Model Estimation via Neural Networks with Statistical Guarantees
Jiahe Lin *
Machine Learning Research, Morgan Stanley
Yikai Zhang
Machine Learning Research, Morgan Stanley
George Michailidis
Department of Statistics and Data Science
University of California, Los Augeles
Reviewed OpenReview: http:://pppaerreeee...eererrrrdd
Abstract
Graphical models are widely used in diverse application domains to model the conditional dependencies amongst collection of random variables. In this paper, we consider settings where the graph structure is covariate-dependent, and investigate a deep neural network- based approach to estimate it. The method allows for flexible functional dependency on the covariate, and fits the data reasonably well in the absence of a Gaussianity assumption. Theoretical results with PAC guarantees are established for the method, under assumptions commonly used in Empirical Risk Minimization framework. The performance of the proposed method is evaluated several synthetic data settings and benchmarked against existing approaches. The method is further illustrated on real datasets involving data from neuroscience and finance, respectively, and produces interpretable results.
1 Introduction
An undirected graphical model captures conditional dependencies amongst collection of random variables and has been widely used in diverse application areas such bioinformatics and the social sciences. It is associated with an undirected graph G = (V, E) with node set V = f1, ,p} and edge set E S V x V that does not contain self-loops, i.e., (j.)) & j E V; x = (21, ... ) is a random real vector indexed by the nodes of V, with probability distribution P.
The connection between the graph G and the probability distribution P(x) comes through the concept of graph factorization and the properties that G exhibits. The probability distribution PP(x) faclorizes with respect to G, if it be written P(x:) 2 ICEε(o (c(#c), where C(G) denotes the set of all cliques of G. PC > 0, VC C(G) potential functions, and Z a normalizing constant. If the underlying graph G in addition exhibits Markov type property (see, e-g-, Lauritzen (1996) and Appendix D.3), then the conditional dependence relationships are encoded in the potential functions pc. In particular, for G satisfying the pairwise Markov property, IP can be written as PP()) [[A,Ee P(j,)) (.)), and the absence of an edge in G implies that the corresponding potential function is zero and consequently the nodes are conditional independent given the remaining nodes. The upshot that the conditional independence relationships for w ~ can be read-off from the graph, factorizes w.r.t. G and the latter is Markov.
For certain multivariate distributions, the form of the joint distribution P(()) enables its decomposition into clique-wise potential functions based on G in fairly straightforward manner; examples include the Gaussian
Equal contribution, ordered alphahetically. Corresponding Author. E-mail: gmichail@ucla.c

---

Published in Transactions on Machine Learning Research (05/2025)
for contimous variables (Lauritzen, 1996, see also Appendix D.1), and the Ising and the Boltzmann machine models (Wainwright et al., 2008) for binary and discrete random variables, respectively. The corresponding potential functions (-(+)) are parameterized by (functions) of the parameters of the underlying distribution (sce example in Appendix D.1). Consequently, a graphical model and the conditional independence relation- ships that G encodes can be obtained by estimating the parameters of the underlying distribution, either through maximum likelihood (Banerjoe et al., 2006; Friedman et al., 2008) or regression procodures using neighborhood selection techniques (Meinshausen & Bühlmann, 2006).
For other classes of graphical models beyond those mentioned above, to select the clique-wise potential functions (c(+) becomes more involved. The difficulty stems from specifying potential functions based on which a compatible joint distribution P(z)) with respect to G can be obtained. One popular approach is to define node conditional distributions through pairwise cliques {j.k}; ie., for any fixed j. P(z (--] x IINEn(()) (),, (x) for k in the neighborhood of j. Such a conditional distribution captures pairwise conditional interactions between nodes, however, this construction docs not directly define proper graphical model for the joint distribution P(ae); unless additional conditions are satisfied (see example in Appendix D.2).
In general, the question of the compatibility of graphical model for the joint distribution (((z) based on a particular specification of all the node conditional distributions in V is rather open, and is generally addressed on a case by case basis (Berti et al., 2014). In this paper, we adopt the convention that models specified through their node-conditional distributions are also referred to as graphical models, which capture node-wise conditional independence relationships.
The presentation thus far focused on a single "static" graphical model, which is suitable for modeling tasks involving homogeneous data sets; i.e., the observed samples æ², 1,.'. ,n are iid from P(æ:). However, in certain applications, the conditional independence relationships may be heterogeneous in that they are modulated by external covariates 2 E RR, An example from biology is that the brain connectivity network- wherein the coordinates of æ correspond to different brain regions-varies varies with age (= E R-). Similarly, gene expression exhibits individual-level variation affected by single nuclcotide polymorphisms (SNPs); see, e.g., Kim et al. (2012). The SNPs are thereby covariates polentially of large dimension that impact the corresponding gene co-expression network. Motivated by such applications, this work introduces a class of covoriate dependent, pairwise interaction graphical models, wherein the log-densities of the node conditional distributions are defined by
J3 = -"
E (5(())o + Sgy fl (1,
(1)
where z a q-dimensional observed covariate, (5x(+) : R9 on z that impacts the strength of the interaction between nodes j and k.²
+ R functional regression coefficient. dependent
1.1 Related Work
We provide a brief overview of approaches aiming to define graphical models for continuous random variables that go beyond the Gaussian case, and also of prior work on covariate dependent graphical models.
Non-Gaussian Graphical models. One line of work focuses on graphical models defined by monotone transformations of the Gaussian case, namely the non-paranormal (Liu et al., 2009; 2012). Their estimation entails first learning in a nonparametric fashion the respective transformations that "Gaussianize" the data, then subsequently leveraging approaches used for the Gaussian graphical model. Another broader line of work focuses on defining graphical models through specification of node conditional distributions. Yang et al. (2014) consider pairwise node conditional distributions from the exponential family, wherein the canonical parameter of the distribution is defined as lincar combination of higher order products of univariate functions of neighbors of each node j in G. This construction leads to family of polential functions that defines a proper graphical model under G for the joint distribution P(x:) (see details in Appendix D.2).
Note that in the case of Caussian distributed error term E, covoriate dependent graphical model is defined in terms of the joint distribution of conditional on 2, i.e., P(]|))) that factorizes over a graph G(z), "parameterized" by the external covariate 2.
2

---

Published in Transactions on Machine Learning Research (05/2025)
More recently, Baptista et al. (2024b) introduce framework wherein the conditional depen- dence/indepeedennce between zj and Tt is characterized by score matrix given by Qzz == 10, Og log p(æ))² where p(z) is the density function of x ~ P(x)) and is further assumed to be con- tinuously differentiable. The authors propose to estimate the density p(z) via a lower triangular transport map, with the components of the map parameterized through linear expansions with Hermite functions as the bases. The same score matrix f is also leveraged Zheng et al. (2023); however, instead of estimating the density p(æ), their approach parameterizes the data score V ((r) using neural networks and estimates it based on an implicit score matching objective (Hyvärinen, 2005), with an additional sparsity-inducing term associated with Ω
Covariate-dependent graphical models. At a high level, covariate-dependent graphical models fall under the umbrella of "varying coefficient models", where model paramcters are assumed to evolve with "dynamic features". Traditionally, the specifie functional dependency on these dynamic features is handled via polynomial or smoothing splines, as well as kernel-local polynomial smoothing (see Fan & Zhang, 2008, and references therein). More recently, it has been rebranded by Al-Shedivat et al. (2020, CEN), where dynamic features are termed "contextual information". Specifically, the authors consider a hierarchical probabilistic framework in which the model parameters depend non-deterministically on the contextual information, and the dependency is parameterized through deep neural networks.
Within the specific context of undirected graphical model estimation where the focus is on conditional independence relationships, the existing literature can be broadly segmented into two categories: (1) the covariate is a categorical variable taking discrete values. Early work by Lafferty et al. (2001) considers a conditional random field (CRF) for discrete æ, and further assuming the underlying graph to be chain graph. For nigh-dimensional settings, with a continuous æ, the problem is cast as efficient estimation of graphical models while taking into account their similaritics across different categorics (Guo ct al., 2011; Cai et al., 2016). (2) For more general dependency on covariates that can be continuous-wal valuud, one line of work focuses on partition- based estimation, namely, segmenting the samples based on the values of the covariates, then performing cstimation within cach partition (Liu et al., 2010). More recently, Gaussian graphical models, in which the graph structure is learned as continuous function of the covariates, are being actively investigated. Specifically, Zhang & Li (2023); Zeng et al. (2024) assume linear functional dependency on the covariate, namely, (z) bjk², z, with b being sparse. The former adopts a regularization- based approach, while the latter a Bayesian one with sparsity inducing prior distributions. Within the Bayesian paradigm, notable contributions also include Ni et al. (2022); Niu et al. (2024) that introduce covariate-dependent prior distributions to enable context-dependent precision matrix estimation. Finally, we remark that graphical model estimation with non-linear dependency on the covariate has been investigated leveraging the idea of kernel smoothing (e.g., Zhou et al., 2010; Kolar et al., 2010). However, these approaches are limited to the case where the covariate is univariate, due to limitations induced by smoothing kernels. Table 1 presents a selective summary of the above-mentioned existing works, with an emphasis on those in which the covariate takes continuous values and corresponding theoretical results provided. We note that the score matrix-based formulation for estimating Man-Gassiian graphical models could potentially be extended to incorporate the dependency on covariates, and such an extension is briefly discussed in Scction 6.
As concluding remark, we note that "contextualized"" estimation has also been adopted in other settings, such as the estimation of directed acyelic graphs (DAGs); see, e.g., Zhou et al. (2022); Thompson et al. (2024) and references therein.
1.2 Our Contributions
This work differs from existing ones in that it investigates a deep neural betwork-based approach for esti- mating covariate-dependent graphical models, and establishes the corresponding statistical guarantecs under the PAC learning framework. Unlike much of the above-mentioned statistical literature that assumes a lin- ear relationship with the covariate, the current work leverages parameterization via neural networks to accommodate more flexible dependency structures. On the theoretical front, finito-sample crror bounds are established under certain regularity conditions, enabling the quantification of the relationship between the size of the neural wetw and the estimation error. Note that although Al-Shedivat et al. (2020) outline
3

---

Published in Transactions on Machine Learning Research (05/2025)
Table 1: Summary of representative works in extant literature covariateedepemdent graphical models.
Zhou at al. (2010)
satup
univariate, Caussian, MLE-based kernel estimator
tre...res.
finito-sample Frohenius and large-deviation bound
Kolar ct al. (2010)
setup
univariate, nodewise regression- based kernel smoothing
theer.
sinite-sample support recovery consistency under sub-Gausanity and minimum edge strength assumption multivariate, Gaussian, tree-based partition & MLE within the region
Liu ct al. (2010)
setup
ree.
guarantee on the excess risk of the estimator
Zhang (2023)
satup
multivariate, Gaussian, nodewise regression w/ linear dependency covariate Anitoo-sample f2 consistency and support rerovery consistency
thacr.res.
al. (2022)
setup
multivariate, Gaussian, likelibood based w/' liuear dependency on covariate no guarantee on the posterior inference of the estimates provided
(Bayesian)
theor.rea rea.
Niu el al. (2024)
setup
multivariale, deparietedepenndell priors & Gauesian liselioood-baee estimationn (wilhin partition)
(Bayesian)
convergence of the posterior distribution in a-divergence metric
This Work
setup
multivariate, nodewise zegression with fiexible dependency on covariate, parameterized by DNN
sinite-sampl fg consistency and edge-wise recovery guarantee with thresholding
Note: for setup, univariate (or multivariate) corresponds to the dimension of the covariste in question, whereas "Gaussian" corresponds to the distributional assumption for in explicitly assumed in the paper.
a generic template for handling contextual information via deep neural networks, the thcoretical results therein primarily focus on (i) quantifying the contribution of the context to predictive performance, under the assumption that the expected predictive accuracy is bounded below by 1- E (Proposition 4), and (ii) analyzing the error bound of the estimated linear coefficient in binary classifications tasks (Theorem 6). Both results are established under the assumption that the underlying true dependency on the covariate is linear.
To this end, the contribution of this work can be summarized as follows.
We consider a nonlinear covariate-depender graphical model, which extends existing ones and allows for flexible functional dependency on the covariate. Empirically, we demonstrate that such an estimation framework effectively learns the dependency of the graph structure on the covariate z; meanwhile, it accommodates complex conditional dependence structures amongst the x and exhibits good performance even when the underlying data deviates from Gaussianity.
On the theoretical front, we establish guarantees on the recovered edges based on a deep neural network estimator. Concretely, a mis-specified setting? is considered wherein both generalization error and an approximation error are present, and their bounds depend on the size of the underlying problem (namely p and q), the size of the neural network and the sample size.
In particular, this work illustrates how existing rosults specifl--ee generalization crror bounds based on local Rademacher complexity analysis (Bartlett et al., 2005) and approximation error bounds of an MLP (Kohler & Langer, (ott)established under generic FRM settings, can be synthesized to obtain error bounds for a specific estimation problem. Such a roadmap is of interest for similar inference tasks.
Note that the roadmap for the technical analysis presented in this paper applics to a broad class of strongly convex, Lipschitz, and uniformly bounded loss functions, by connecting the infinity norm of the function approximation in Kohler & Langer (2021) to the excess error bound in Bartlett et al. (2005) through a strong convexity- argument. Hence, it goes beyond the result established for the mean squared error loss function studied in Kohler & Langer (2021).
2 Dased Govariate-Dependde Graphical Model Estimation
We elaborate on the model under consideration in (1) and provide its formal definition next. We consider a system of p variables == Z2, , ,zp)' taking values in X C R2, whose conditional dependence structure depends on some g-dimensional "external" covariate E Z C Rq. Let [p] {1, (1....,p)) ,p} denote the index
Here mis-specified in the ERM context, where the target hypothesis does not necessarily live in the hypothesis class in which empirical minimization conducted.
4

---

Published in Transactions on Machine Learning Research (05/2025)
set of the variables. To model such conditional dependence structures, consider nodewise regression formulation as the working model, given in the form of
Jj - (8,(2),
+ Σ B,(())n t e lv].
(2)
E(E,) = and has finite moments; (-x 2 + R is a function that takes & as the input and outputs the regression coefficients for v 's's with 2; being the response. For notation simplicity, we let 8)(2) 2 RP-1 whose output coordinates correspond
k -, -...... .P, + j, and æ-j = (ø1.-... We note the analogy to the node-wise regression formulation for the Gaussian graphical model (Meinshausen & Bühlmann, 2006) (see also Appendix D.1 for brief discussion), where in the absence of the external covariate the regression cocfficient 35x "degenerates" to a scalar. With the dependency on 2, the graphical model is solely captured in the (z)'s, and one can retrieve the estimated graph at the sample level, based on the corresponding & value for the sample of interest. In particular, when æ is jointly Gaussian, the M(Eie, Bis (2)
node-conditional distribution is given by ((rz
We consider parameterizing each (-x(-) by some neural network, e.g., a multi-layer perceptron (MLP). Let o denote the parameters of the neural networks in question collectively. Given independently identically distributed (ii.d.)) samples [(x', >²) 0 can be obtained by minimizing the empirical loss; e.g., with a mean squared error loss, it is given by
arg min
| x' 112 where 8) (B,(==))),
(3)
the dependency on the neural network parameter 0 is made explicit. The training pipeline is outlined in Exhibit
Exhibit 1: DNN-based Covariate-deppnnentt Graphical Model (DNN-CGM) Learning Pipeline
Input: Input i.i.d. samples (()²==
1. [forward pass] for fixed 0, calculate ;) for all j = 1, .................. p and i = 1, .............................. according to (3);
2. [loss] calculate the loss based on all samples, by evaluating the distance between x' and a';
3. [backward pass] update 0 using the gradients calculated based on the loss (back- prop-proppg:t:on)::
Output:
Estimated 1 and (B,(21, 0)
Remark (On the implementation of (.x networks). Instead of creating separate neural works for each ())(+) which requires p(p-1) networks, with each taking z as the input and generating a scalar output one can alternatively create neural network(s) with appropriately sized hidden /output layers that "backbones" are shared across and fewer networks are needed. This leads to fewer total number of parameters needed and empirically an easier optimization problem, which often yields superior performance. As a concrete example where a single neural network is used, it takes z as the input and outputs a p(p 1) dimensional vector, with each coordinate corresponding to the value of (z) j. k E [p]; k # j: such a design corresponds to the case where 3z)(2) operationalizel as (2) = (u))). Specifically, let h be the dimension of the last hidden layer; u(-) : IR" + RR is shared across all 8yx's and corresponds to stacked neural network layers up until the last hidden layer; (-)(-) R² + R is linear layer that maps neurons from the last hidden layer to the output unit indexed by jk, and it differs across the jk's. In this case, 0 = (0., (Or.jk: k E [p]: k j]). Remark 2. The formulation under consideration not restricted to z being continuous. In the case where any coordinate(s) of is categorical, one can first process the categories via an embedding layer, then concatenate them with the numerical ones and proceed.
Some considerations on the formulation are discussed next. First note that for a node-condiiional formulation where the conditional distribution of cach node is some function of its neighborhood set (Meinshausen & Bühhmann, 2006), the model in its most general form can be written as zÿ = f=(z=),zz)) E jp). Such a form can be too general, and potentially incur difficulty in actually identifying the neighborhood set based
5

---

Published in Transactions on Machine Learning Research (05/2025)
data?. In the same spirit in Büihlmann et al. (2014), by considering an additive form, the model is restricted to z-j = Exyz [s((xs) + ej; in particular, lzx(zz) 0 for k & ne(j); ne()) S [p] \ (5) is the neighborhood of node j. With the presence of the covariate = that impacts the neighborhood structure, a natural extension is to augment the input space of fzx s,, namely,
B3 T= = Mey
Ent ()((2, tey e [p).
In this paper, we impose additional structure on the specification of fyk and let f5((2, = Bz((z)=z. The implication is two-fold:
1. z and Ih are assumed separable.In particular, we consider a multiplicative form, namely, (z, xR) = (z)())); that behaves like a "locally linear" model on Y5k (xx), and 55x(2) can be viewed as a generalized regression coefficient. The edge weight can be summarized by B)x(2), which can be This pre-
interpreted as the strength of relevance; in addition, 35x(2) = = fzz(2, mx) serves the interpretability of the model in that the neighborhood set can be solely inferred based on {=)=(<)))=1,->>> pkksj, a primary quantity of interest that can be directly estimated and prob- abilistic guarantee can be established accordingly. Note that for fzx(z, that allows for general interactions between 2 and The extracting the underlying weighted graph edge that characterizes the strength of the dependency requires past-hoc operations. For example, this can be achieved by considering quantities of the form (=z,).22
A similar separable structure is employed in Alvarez Melis & Jaakkola. (2018); Marcinkevičs & Vogt (2021), where, despite a slightly different model setup that involves only æ, the function is decomposed as ((z)h(x),, where h(x) serves as a representation of x and 0(x) its corresponding relevance; those models are interpretable and "self-explaining""
2. The functions (rx)ss are further simplified to (Tk())) = dr. From a regression analysis perspec- live, Iy and its transformations can be interpreted as features. We use "raw" features directly, partly due to the fact for the graphical model setting under consideration, they correspond to node values and thus are inherently meaningful. This contrasts with, for example, image pixels that lack intrinsic intepretability and thus benefit from transformations into higher level features (Alvarez Melis & Jaakkola, 2018). One can potentially consider a more complex (5x()), e.g., pa- rameterize it by some neural network, or through basis expansion (e.g., Qiao et al., 2019), namely, 75x (zx(2) == Ela (x)) with the bases Φe(-) fixed apriori and ajk,e learnable parameters. Empiri- cally, we observe that the current specification is essentially as powerful as its enriched counterparts (modulo tuning), in terms of identifying the neighborhood of each node.
Remark 3. Up to this point, we do not impose any sparsity assumptions on the underlying graph and the estimation procedure does not involve any regularization terms. In practice, to obtain a sparsified graph for interpretability purposes, we consider a thresholding procedure, where the small entries can be effectively treated as zeros. The thresholding level is selected by examining the edge magnitude histogram, and chosen in the region where a "gap" is present Finally, given the undirccted nature of the graph, we adopt the "AND" principle, namely, we consider an edge as present when both Bye and (ay are nonzero after hard-thresholding.
3 Theoretical Results
This section provides theoretical guarantees for the model posited in (2), assuming that it corresponds to the ground-truth data generating process (DGP). Specifically, let (*)-)) denote the target model (i.c., the true DGP one); given the availability of data ((x', , generated according to (2), we are interested in the
"Wtt a general fj where interactions amongst the zs's are permitted, for æ; IL Fg x- to hold, z7 needs to be absent from all terms, which is difficult to operationalize either through constraints in post-hoc fashion.
SSe also the setting considered in Corollary 2. Note that in the case where the underlying true graph has edges that can be separated into weak and strong ones (with exact sparsity being special casc), one would expect gap in the histogram, aside from estimation/ cstimation/approcinaate /aprroximation errors. Empirically, one can first "normalize" the graph so that the maximum entry does not exceed in magnitude; the thresholding level is then scaled accordingly. This normalization step can potentially enhance the interpretability of the the thresholding level, without fundamentally altering the sparsification outcome.
6

---

Published in Transactions on Machine Learning Research (05/2025)
error bound of (()) under the Empirical Risk Minimization (ERM) framework, with ()x parameterized by deep neural networks. We establish finite-sample bounds for the generalization error using Tocal Rademacher complexity tools (Bartlett et al., 2005), and the approzimation error bound in terms of neural network hyper- parameters (namely, depth and number of neurons) leveraging the techniques in Kohler & Langer (2021), under certain regularity conditions. The two combined provides insight into consistency type of results for B()), relative to the ground truth ()).) In the ensuing technical developments, we further assume that the random vectors 2, Z are bounded, namely E [-C, CPP and CJC respectively, Without loss of generality, we set C 1, which only affects the scale in the analysis.
3.1 Preliminaries and Road Map
Definitions. Let (), P) denote probability space and Ey[-] denote the expectation with respect to P over random variable y V) Let H denote hypothesis class; i.e., class of measurable functions h y + R. The infinity norm between two functions h, h is defined as j/h-hlloc: := supyey Ih(y) - -h(y), while for vector functions h = (h.. ,hp), the corresponding vectorized infinity norm as ||h hl x == max 1| hy bllo- For two functions, h, g. the notation h Z means that h 2 cg for some universal constant c, and S is analogously defined. Further, g = h, if cg S h S ēg for some universal constants ē and Finally, for a real-valued vector y, its weighted norm with respect to a matrix A is given by ||yl == Ay).
For the problem at hand, let y := z); Sn == {(x', (<=',²1)) >))) C (X x Z)" denote a random sample of iid data points. We consider a loss function b(a;b) : R + and H. the hypothesis class of risk functions of the form h(y) == (y) == 1 (E kzj fix(2) (x)),,,)))) i.c., the risk is a composite function of the loss and regression coefficients. Then, the population and empirical versions of the risk are respectively defined as follows, with [n] = f1. .n]:
E
(18,(t), (1))),
R(8) = Fa.s
(4)
(B) = Σ [m(m) (10(())))
To perform ERM, we assume that the regression cocfficient functions (-x(-) 2 + R belong to some hypothesis class Fo; and denote Fpx(p-1) as the joint of x (p - 1) hypotheses, with each individual hypothesis B-x(-) E F,j,k lp), k + j. We further define collection of optimal hypotheses 87
gopt's as
53 == arg s eF 18.e Bjulloe, Vj. [p]: #
and the approximation error can be defined as Expprox (F) = maxje(1," vi) l18 8illoc. Finally, the ERM estimator the regression coefficient functions is given by
D == arg miggiinguux (3).
(5)
Assumptions. Before presenting the results, we posit the assumptions used in the theoretical analysis. Together with the boundedness assumption on (x,z) mentioned above, these assumptions lead to a O(1/n) fast convergence rates for the risk excess error, and hence a smaller sample size requirement. Further, the analysis enables us to balance the trade-off between approximation error and generalization error based on the size of the neural network, the sample size and the dimension of (x)
A Assumption (On the loss function). The loss function E(a; R R [0, w.o.lg. we set M 16. We assume E(a, b) is:
I--iipschila in a: (l(a1:b) f(a2; b)l < L|a1 - 02l, where L > O denotes the Lipschila constant.
a-strongly convex in a: E(a1;b) E(oz; b) 2 (((aa) a=a (a1 a2) + (a1 o2)² where 0 and (()-) denotes the derivative of the function l.
5 is typically selected by practitioners, such neural networks with fixed number of layers and neurons.
For any loss function ((a;6) Rx R + (0, M]. suffices to consider rescaled version 1. -((a; b) to satisfy such assumption.
7

---

Published in Transactions on Machine Learning Research (05/2025)
Remark 4. Assumption 1 is widely used in the statistical learning literature, especially when aiming to establish last rates for the generalization error bound (see, e.g, Klochkov & Zhivotovskiy, 2021).
Assumption 2 (On the optimality of (.x)- The optimal regression coefficient function 3J.5 [p] satisfies the following condition:
Le,z Th
= 0.
Remark 5. Note that the combination of Assumptions 1 and 2 implies that 5. E [p], k + minimizes the population risk in (4) due to the strong convexity of the loss function l. We also demonstrate that the mean square error loss function satisfies Assumptions 1 and 2 in Appendix A.2.
Assumption 3 (Bounded pseudo-dimension of coefficient functions Bj. E [p] has finite hypothesis class F). The hypothesis class F of the regression (Pollard, 1990); i.e., dp(F) < oo (see formal definition of diauddd-mmensionn in Definition 3 in Appendix A.1).
This assumption is widely adopted in the statistical learning literature and encompasses broad range of functional classes, including linear and polynomial functions (Anthony & Bartlett, 1999), as well as various families of neural networks (Bartlett & Maass, 2003; Bartlett et al., 2019; Khavari & Rabusseau, 2021).
The primary result established for the model postulated in (1) is the statistical consistency for the estimate 3 (Corollary 1). This is achieved through two key steps: (i) derive bound on the excess error, defined as ((8))) as function of (3; Theorem 1; and (ii) leveraging the strong convexity of the loss function (mean squared error loss) and the result in Theorem 1, derive the error bound for under some projection norm. The bound in Theorem 1 has two components that respectively capture the generalization error and the approximation error, with the latter arising from the mis-specified setting assumption that allows the truc 9 to live outside of the working hypothosis class F. To derive bound for the generalization error, we first show that the hypothesis class that encompasses composite functions quantifying the "delta in risk" between & and the true 9 satsifisies Bernstein's Condition (Lemma. 1). Next, we select an appropriate subroot function r(-) tailored to the problem, enabling us to apply the proof strategy and results from Bartlett et al. (2005) (specifically, Theorem 3.3, Corollary 2.2, and Corollary 3.7). A critical component of the proof involves bounding the covering number of the function class 35x (z) under consideration.
With the working hypothesis class encompassing class of MLPs, Theorem 2 characterizes the approximation error for such architectures, employing tools from Kohler & Langer (2021). Note that the pseudo-dimension which appears in the generalization error bound- can also be associated with the size of the neural network. Consequently, the two error components for the MLP architecture under consideration can be "balanced", in that one can minimize the excess error bound by choosing the quantity that links the two appropriatcly.
3.2
Main Results
The first result leverages Local Rademacher Complexity and empirical processes tools to bound the gener- alization crror:
O dP(FF
R(3) (R(BR()) (
excess eITOr
empirical excess error generalization error
where the O(-) notation ignores all problem dependent constants and only depends and dp(F). The empirical excess error can be further bounded by an approximation error term controlled by the capacity of the hesis class (see equation (15) in the proof of Theorem 1 in Appendix A.3)
Theorem 1. Under Assumptions 1-3, the following bound for the excess error as a function of the empirical
risk minimizer 99 (sce (5)) based on S, samples holds with probability at least 1 5, for any d >>> 0:
Bdd(F)
R(Ø) - R(B)) L- Eapprrxx (F)
log(nLL)) log
(6)
8
approximation error
generalization etvor
8

---

Published in Transactions on Machine Learning Research (05/2025)
Remark 6. The p² term in (6) shows up due to estimating p x (p- 1) functions using neural networks, which naturally increases the entropy number of F with respect to S, in Dudley's 's entropy integral (Dudley, 2016) by a factor of p². This indicates that the size of the graphical model can grow at a p o(vn) rate for the excess error to vanish. A larger p can be accommodated, cither by incorporating sparsity assumptions (Meinshausen Bühlmann, 2006; Wainwright, 2009) or through refined analysis as multi-lask learning (Lounici et al., 2009), which however is outside the scope of the current paper.
The following corollary establishes the consistency of the empirical risk minimizer, leveraging the strong convexity the loss function.
Corollary 1 (Consistency of 3). Let A₃(2) =- Le (=- |=|. Under the Assumptions of Theorem 1, the following inequality holds with probability at least - 8 for any > 0:
P (1,(2)) -3(=))2)))))) <
LFd
p² log(nLp) log (1)
L Eapprox (F)
a²n
Corollary 1 quantifies the quality of the recovery of the regression coefficient functions 33 that capture the strength of the conditional dependence relationships between the variables.
Note that to establish Theorem the assumption of realizability i.e,, (+) E is not imposed. Con- sequently, the empirical risk of the ERM estimator is not guaranteed to be smaller than that of B i.e., (B*) & Rn is not guaranteed. In particular, this can happen when the hypothesis class does not have adequate capacity, and hence the approximation crror term is nonzero and appears in (6). However, Theorem 1 has not. characterized the approximation error terms of the sample size and dimension p of the graphical model. This is established in Theorem 2, which also aims to balance the approximation error and generalization error terms.
To that end, in addition to Assumptions 1-3 used to derive the excessive risk bound, need the following additional assumption to apply tools from Kohler & Langer (2021) to quantify the approximation error based on the number of layers and the total number of neurons used by the feed-forward neural networks to approximate the regression coefficients 8yk-
Assumption ((m, C)-smoothness of (5 We assume 87. is (m, C)-smooth for all j, E [p), + (see the formal Definition 4 in Appendiz A.1).
Theorem 2. Let F be family of fully connected neural networks with ReLU activation functions, hav- ing number of layers H = g- and number of neurons = (2e)9 q g². Under the Assumptions of Theorem 1, together with Assumption 4 with S q, by setting
[.mm'd log2 (mIp) log(p) log(1 (1)) log(1/)) m+ {= {=(
an
the following holds with probability at least - 8:
R(A) R(()) $ T.g and Fg 1 E 3=1 110,(2) -03(()))aaaa))
S LE
(8)
Remark 7. The rates of the bounds in (8) are of the order n mmmm(m+ whereas the rate in Kohler & Langer (2021) based on a quadratic loss is n 2mm (2m-1), This discrepancy is due to the fact that when analyzing general I-TTisshhiz functions, the excess error bound depends a linear fashion on the approximation error of 8(z). For the mean squared error loss function, it depends in a quadratic fashion and hence the same rate as in Kohler & Langer (2021) can be obtained, as shown in Appendix A.6.
Extensions to edge-wiss wise recovery. Recall that in Corollary an t2-type consistency result is established for 3 based on a projection norm. If we consider a setting where edges can be segmented into strong and weak ones, then the result in Corollary 1 can be translated to an "edge-wise" type guarantee, Formally, we define the set of strong odges as E (z) == {G. k)] 185x Z 2) and of weak ones as E=)(z)= == {(j,A) (0)%( s 3). where 3 and correspond to their respective catiiimmmmaiimmni magnitude. In addition, we assume that
9

---

Published in Transactions on Machine Learning Research (05/2025)
there exists uniform lower bound Φ 0 on the margin between the maximum and minimum magnitudes; namely, _ B 2 > 0. The edge-wise result. in Corollary 2 focuses establishing guarantee that
(=)3 2 matches ([]] 2 {} in expectation, for some threshold T = n/3 + (1 - 1) " > 0; i.e.,
is convex combination of 3 and B. Assumption 5. = ((==)==[[=_a=_])_
yI holds wniformly for all 2 E Z with y > 0.
This assumplion ensures that the minimum eigenvalue of A₃(=) is bounded away from zero, uniformly for all This assumption allows to to first establish an analogous result as in (7), with the distance between (2) - (=) measured in the Euclidean norm || instead of the projection norm Az):); subsequently, the Euclidean norm result is further translated to that on individual edges.
Corollary 2. Let T == v1 + (1 - m) with 7) > O. Under the Assumplions of Corollary 1 and Assumption 5, the following inequality holds with probability at least 1:
Σ Σ (1 2 20} j=1 And
[.d
p² log(nLp) log (1)
a² min(1r (1 -) 13)333)3)
(9)
L. Eapproxx (F) ary mymin/17 (1 n) (3 3)2
Remark 8. Some comments on the result in Corollary 2 are provided next.
Similarly to the result in Corollary 1, the bound in (9) consists of two terms that correspond, respectively, to the generalization error and the approximation error bounds. With a constant margin 3- the generalization error bound still vanishes, provided that the network size and the sample size n grow at certain rates. However, the approximation error bound does not vanish; in particular, when the margin is small, the corresponding term can be large, resulting in a fairly loose bound. This aligns with intuition that when the strong and weak edges are hard to "distinguish" the edge-wise recovery via thresholding becomes difficult.
The established bound is in expectation and can be interpreted as follows: on average, the strong edges can be identified via thresholding, provided that the threshold is selected within certain range. This differs from the graph recovery results established for sparse high-dimensionall Gaussian graphical models based on an & penalty (e.g., Ravikumar et al., 2009), largely due to the fact that we are under PAC learning framework and hence there is some discrepancy embedded in the setting in question-thar we consider a mis-specified setting without assuming that we know the functional class of the underlying true data generaling process, and thereby the tools adopted.
An important implication of this result is that, in practice, when there clear separation in mag- nitude between strong and weak edges, one can effectively obtaín a sparse graph by thresholding. This approach practically treats the weak edges as zero, thereby enhancing interpretability.
Remark (Refinement of Corollary 2). By incorporating the analysis on the functional approximation error, the above corollary can be refined as follows, under the same setting as considered in Theorem 2:
E.
{f 2 * (8A(())
2
PLE 21] & min 7)² B)²
4 Synthetic Data Experiments
The performance of the proposed ncural network-basod mothod is asscssed through a scries of experiments on synthetic data sets. Both Gaussian and non-Gaute non-Gausiian settings are considered, with samples (indexed by i) generated according to one of the following three mechanisms:
Gaussian: x' (O, (0))-)).
Non- con-parnnormal (NPN, Liu et al. (2009)): x' = f- ((=); x² ~ ~NO0 (e+)-)) with f being a monotone and differentiable function, and its inverse f-1 is applied to o² in coordinate-wise fashion.
10

---

Published in Transactions on Machine Learning Research (05/2025)
Directed Acyclic Graph (DAG): for each coordinate z = Exepatr) f3((t)) c,, where pa(j) is the porent set of node j. In other words, x' is generated according to a structural equation model (SEM). Note that the SEM serves solely as a data generation mechanism to introduce polentially nonlinear dependencies (e.g., by parameterizing fas's through basis functions) amongst the nodes in flexible way. However, the primary quantity of interest is still the conditional independence structure captured in the corresponding undirected moralized graph, as described next.
For the Gaussian and the non paranormal cases, O' is determined by z'. For the DAG case, we start with a binary A determined by z' that corresponds to the skeleton of the DAG (denoted by DAG); the parent set of node j satisfies pa(i)= = {k: : Ais + 0}. The undirected conditional independence graph of interest that is compatible with the data can then be obtained by moralizing the DAG (Cowell et al., 1999, Chapter 3.2.1), namely, G aa))) whose graph structure is encoded in O It can be viewed as the counterpart to those in the Gaussian haassian/non-Paa moon-paranormll case in terms of capturing conditional independence relationships". For all three cases, Ot is the parameter of interest, and we are interested in how well its skelelon can be recovered by the proposed method. The results are benchmarked against various competitors, including RegGMM (Zhang & Li, 2023), glasso (Friedman et al., 2008) and nodewise Lasso (Meinshausen & Bühlmann, 2006). RegGMM accounts for covariate dependence while assuming linearity, whereas the other two are graphical model estimation methods that solely consume æ as the input.
Settings. The description of the various simulation settings is outlined in Table 2. For settings where data are generated according to the Gaussian or NPN models, we first generate "candidate' precision matrices denoted by Vt's, where V either has a single band steps from the main diagonal (for settings G1/N1) or is block diagonal with nonzero entries on the (-th block (for settings G2/N2). Next, each O' is a conver combination which automatically ensures the positive definiteness of O of the candidate precision matrices, and the mixing depends the value of the covariate x². For settings D1 and D2 where data are generated according to DAG, we first generate binary matrices B1 and B2; both B1 and B₂ correspond to the skeleton of some trees, with nodes having 1 to 3 children (randomly determined as the tree grows). The covariate z dictates the skeleton of DAG, i.e., A", and also governs fix either directly through the coelficients on the zg's 's (linear case), or through a multiplier that impacts the coefficients of the basis functions (non- lincar case). For all scttings, the samples effoctively fall under different "clusters" according to their z"ss. Samples within the same cluster have identical skeletons, albeit the magnitude of the entries may differ depending on the exact value of the z's.
The degree of lincarity in z varies by setting. For scttings G1 and N1, the dependency on z E R² is linear in its 1st coordinate which governs the mixing percentage, and the 2nd coordinate ellectively dictates the cluster membership. This is similar for settings DI and D2, which however differ in that the dependency on the 1st coordinate being quadratic. For settings G2 and N2, the dimension of 2" is set at 10; to induce non-linearit, we use radial basis function (RBF) network &p that first transforms x into a scalar, followed by a sigmoid transformation, namely gi == sigmoid((<(22)) where p(z) == Ete ag cxp(- - cxll²); the sigmoid transformation ensures that # E (0, 1), which then dictates the mixing percentage. Finally, note that amongst these settings, the dependency on the ay 's are linear in settings GI, G2 and D1. Additional details and pictorial illustrations are deferred to Appendix B.
For all settings, we train the model with B(-) parameterized with an MLP on 10,000 samples, and evaluate it on a test set of size 1000; the estimated graphs (off-diagonals only) for these test samples (indexed by i) are extracted as { (=3))(2')) (,)))) For benchmarking methods that are not covariate-dependen. i.c., glasso and nodewise Lasso, their sample-level estimates are identical, namely θ² = θ for all i. Instead of directly running the mothod on the full set of training samples, we further partition them based on their
II the special case where is generated according to linear SEM, namely x = Az + e, e == cov(2) satisfies θ-1 = -T where Q = cuv(e); see, e.g, Loh & Bühlmann (2014). In the absence of lincarity. the exact magnitude of the entries in @ becomes difficult to infer; however, one can still infer its skeleton from that of the DAG by the definition of moralization.
(I - AT) T
TT( -
We use stacked linear layers with ReLU activations, and do not leverage any RBF-style layers or sigmaid ones which are used the data generating pracess.
11

---

Published in Transactions on Machine Learning Research (05/2025)
Table 2: Overview of simulation settings.
MECHANISM
DESCRIPTION
41 + (1 - 2)
if c (0,1/3)
Gaussian
p 50.g = 2:
Unif(0, 1); @²
#1 V2 (1 +
if c (1/3,2/3)
V1 + (1-
p 90,9 q 10; x² ~ standard multivariate Gaussian; 8 sigmoid((<=*)); sp()) is an RBF network;
x0z) = Σ 04 exp(p oxll?); L= 10, or ~ Unif(-10, 10), 8r ~ Uni or ~ [Unif[- Unif(-1111)" 1),4
G2
Gaussian
/ + v3
if if 0.9 or 2 < 0.1
=
// 0.5 1 (0.5 Va otherwise
N1
NPN
is's are generated identically to those in GI; x' = g(²)) where g(x) = x + sin(x).
N2
NPN
i are generated identically to those G2; x² = g((*) where g(x) sign(()).
B1
if o C (0
p = 50, 2; Unif( 1): A'=
B2
= 6- (-2.0);
D1
DAG
otherwise
(a(a) = #x: A == 1( #0).
g. zts and At are set identically to those in D1; f56(2) == 1
v1 (2) a)a2)z(a) + 33 (xx(2) 10m (-'s are Hermmte functions (Olver et al., 2010); a ajk. = (A%, -.a.m). ~ Unif(0.1, 0.5).
D2
DAG
cluster membership, and conduct separate estimations using only within-cluster samples. 9 Note that such partition would not be feasible in real world settings as the cluster membership would not be known apriori. All experiments are repeated over 5 data replicates.
Performance evaluation. We use AUROC and AUPRC as metrics, which adequately capture how well the methods estimate entries with strong signals versus weak ones"", and there is no need to apply hard- thresholding to estimates to obtain these two metries. For a single experiment based on (any) one data replicate, the metrics are initially obtained at the sample-lovel by comparing O' against 0i for all samples under ";nsideratiou to oblain the metrics of interest corresponding to a single experiment, we average those over the samples. For glasso /nodwwise Lasso, experiments are conducted over a sequence of penalty parameters, and the highest metric values are reported. Table 3 reports these metries, after averaging across the data replicates (experiments), with the corresponding standard deviation reported in parentheses.
Table Evaluation for the proposed and benchmarking methods, averaged over 5 data replicates
DNN-based CGM
RegGMM
glassso - est. by cluster
nodewise Lasso est. by cluster
AUROC
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
AUPRC
GT
0.59 (0.000)
0.98 (0.001)
0.36 (0.002)
0.79 (0.004)
0.97 (0.000)
0.59 (0.000)
0.33 (0.000)
0.98 (0.002)
G2
0.99 (0.003)
0.92 (0.024)
0.86 (0.017)
(0.020)
0.97 (0.000)
0.45 (0.031)
0.96 (0.007)
0.55 (0.141)
N1
0.99 (0.000)
0.98 (0.001)
0.99 (0.001)
0.89 (0.001)
0.98 (0.000)
0.59 (0.001)
0.99 (0.001)
0.98 (0.003)
N2
0.99
(0.004)
0.90 (0.028)
0.82 (0.015)
0.37 (0.023)
0.96 (0.000)
0.42 (0.031)
0.93 (0.016)
0.37 (0.16-1)
0.95 (0.004)
0.76 (0.026)
0.94 (0.005)
0.74 (0.026)
0.92 (0.005)
0.46 (0.008)
0.94 (0.005)
0.85 (0.032)
D2
0.92 (0.008)
0.66 (0.022)
0.87 (0.013)
0.60 (0.017)
0.89 (0.009)
0.10 (0.017)
0.87 (0.015)
0.50 (0.113)
The main observations are: (1) the proposed method exhibits superior performance than existing bench- marks, and particularly so in AUPRC. As matter of fact, all methods exhibit reasonable performance in AUROC, even for /nodwiiee Lasso that cannot perform sample-specific estimation. Notably, when the dependency on z is linear within each cluster (scttings G1 and N1), nodewise Lasso performs almost perfectly when separate estimation is conducted within each cluster and thus on samples whose underlying
PReall that that the true graphs corresponding to samples within the same cluster have identical skeletons, although the magnitude may still vary. glass/noddwwwr Lasso are still "mis-specified", ; however to a lesser extent when compared against the case where this "cluster membership" information is ignored.
10The estimated 9 does not correspond edge probabilities, but rather, it corresponds to edge weights. AUROC/AUPRC can effectively summarize how well the graph skeleton is being captured after the entries thresholded over range of thresholding levels.
11 For methods that can produce sample-ppecific istimates, the evaluation is done the test set; for those that can only produce a single estimate based on all samples, the evaluation is done on the graph estimated based on the training samples directly.
12

---

Published in Transactions on Machine Learning Research (05/2025)
true graphs have the same skeleton yet different magnitudes, despite that the model assumes that all samples have identical graphs. (2) In the case where the true data generating process is mildly non-linear in zy in'siin particular, the non- non-linearity is induced via monotonic transformations (i.e., NI and N2), the proposed formu- lation in (2) as a working modcl-by only considering linear dependency on the zy. 's is still able to recover well the skeleton. (3) In the case where the true DGP becomes highly non-linear (e.g., D2), the performance of the method deteriorates. This expected in that the moral graph can be inherently difficult recover, partly due to the existence of edges that are induced by "married" nodes and have small magnitude; the non-linearity in I&'s 's further add to the complexity of edge recovery. The former is corroborated by D1, where although the dependency on Z& 's is linear, the skeleton recovery is still inferior. Additional tables and remarks are deferred to Appendix B for further benchmarking and illustration, including the skeleton recovery (0/1) performance DNN-CGM at different thresholding levels (Table 4), the performance for the DAG setting when estimates are evaluated against the "pseudo" moralized graph without the edges from the married nodes (Table 5), and the performance of glasso/ nodowisc Lasso without sample partition (Table 6).
5 Real Data Experiments
To demonstrate how the proposed method performs in real world scttings, consider two applications from the domains of neuroscience and finance (deferred to Appendix C), and assess the network structures recovered.
We consider a dataset from the Human Connectome Project analyzed in Lee et al. (2023) comprising of resting-state IMRT scans for 549 subjects. The scan of any subject in the form of a spatial-tempon matrix, with the rows corresponding to the "snapshot" value of 268 brain regions and the columns being the temporal obscrvations of the region, totaling 1200 time points. In addition, cach subject is associated with Penn Progressive Matrices score, surrogate of fluid intelligence. The quantity of interest is the network of brain regions (p- 268) function of the intelligence score (q 1). The dataset has been pre-processed with global signal regression, where shared variance between the global signal and the time course of each individual voxel is removed through linear regression (Murphy & Fox, 2017: Greene et al., 2018). In our experiment, we flatten the time dimension the IMRT scans and ignore the temporal nature of these scans; therefore, the dataset contains 1200 549 660k effective samples. However, note that since the score is the subject level, the estimation procedure yields single unique network for each subject
quartle-qt deretyy (1.dr=11.41
Startllwugt denutyy (utweigntest-221
Lorrz
Unoond
(a) Network for subjects whose scores in the 1st quantile (b) Network for subjects whose scores are in the 4th quantile
Figure 1: Estimated networks as represented by -3)) for subjects having low (left) and high scores (right), after thresholding at 0.05. Red cells indicate positive partial correlations and blue cells indicate negative ones.
12This due the fact that model specified function the score and thus effective samples of the subject will always have their sample-specific estimate graphs being identical.
13

---

Published in Transactions on Machine Learning Research (05/2025)
Results. We group the subjects based on the quantiles that their respective scores fall into, and obtain the averaged graph over subjects that are within the same quantile band (labeled as ql, q2, q3 and q4). Further, for visualization purposes, we partition the nodes into 6 sub-networks, corresponding to default mode, frontoparietal, medial frontal, motor, subcortical cerebellar and visual, as in Lee et al. (2023). Figure shows the comparison between the medial frontal sub-networks of subjects in ql and q4, respectively, after the estimated graphs are sparsified (see Remark 3) and thresholded at 0.05 (on a normalized scale). Note that each node is associated with label's, in the plot, nodes with the same label are grouped together.
The difference in the network density is rather pronounced; in particular, those having a higher score exhibit significantly more connected network. Similar patterns are also observed in other sub-networks; in particular, the differential between qI and q4 networks frontoparietal, motor and visual is of comparable scale to that in medial frontal, while that of default mode and subcortiical-cerebbllla. is less pronounced. This is largely concordant with observations from existing litcrature (Song et al., 2008; Ohtani al., 2014) and is also corroborated in validation dataset encompassing 828 subjects; the estimated graphs therein exhibit the same pattern.
6 Discussion
A nonlinear covariate-dependent graphical model based on nodo-conditional formulation is investigated. The functional dependency on the covariate z is parameterized by neural networks, and hence the model can capture this dependency in flexible way. Theoretical guarantees are provided under the PAC learning framework, wherein both the generalization error and the approximation error are taken into account.
At the methodological level, allernative approaches similar to the formulation in Baptista et al. (2024b), can be potentially adopted, by considering a modified score matrix that takes into account the external covariate, namely Ωy* [](=]= log [[x]==]]² Consequently, estimation can proceed with the use of a lower triangular transport map to estimate the conditional density p(z x() (Baptista et al., 2024a), score matching to approximate the conditional score Va p(a)=)) (Dasgupta et al., 2023). However, it is worth noting that the consistency result established in Baptista et al. (2024b) is under an asymptotic regime, where the number of nodes is assumed fixed and the sample size n grows. Their analysis relies on a Taylor expansion, followed by steps that invoke the delta method and the continuous mapping theorem, under some additional assumptions In contrast, the thcoretical result in this paper establishes a finito-sample error bound for the estimator under regime where p can grow slowly with n. Hence, the analysis requires different set of technical tools, even in the absence of the covariate z. By adotping node-conditional formulation, the learning task reduces to an ERM problem, enabling us to leverage existing results in the supervised learning literature.
Finally, the focus of the paper is on graphical models for continuous variables æ, but with appropriate modi- fications the proposed framework can be extended to discrete random variables e.g., a covariate dependent Ising and/or Potts model (Wainwright al., (008)-with a cross-entropy loss function. Theoretical results can polentially be established using a similar set of arguments, leveraging results in ERM for the case of classification.
Acknowledgements
The authors thank the Action Editor and three anonymous reviewers for their careful review of the work, and their constructive comments and suggestions.
George Michailidis was supported in part by NSF grants DMS-2348640 and DMS-2319593.
13 mongst all 268 nodes, there are 20 distinct labels
1d Two key assumptions include: (i) the parameterization based on bases expansions of the transport map sufficiently rich to cover the target density, and (II) the estimated before thresholding an exact expectation, namely, Q,a == El0,ö, log p())) rather than sample-level estimates.
14

---

Published in Transactions on Machine Learning Research (05/2025)
References
Al-Shedivat, M., Dubey, A., and Xing, E. Contextual explanation networks. Journal of Machine Learning Research, 21 (1914)):144 44, 2020. URL http://jmlr arg/paaprrsv21/11--022 html.
Alvarez Melis, D. and Jaakkola, T. Towards robust. interprelability with self-xxlaining neural nelworks. Advances in Newral Information Processing Systems, 31, 2018.
Anthouy, M. and Bartlett, P. L. Neural Network Learning: Theoretical Foundations, volume 9. Cambridge University Press, 1999.
Banerjee, O., Ghaoui, L. E., d'Aspremont, A., and Natsoulis, G. Convex optimization techniques for fit- ting sparse Gaussian graphical models. In Proceedings of the 23rd International Conference on Machine Learning (ICML), pp. 89-96, 2006.
Baptista, R.,, Marzouk, Y., and Zahm, O. On the representation and learning of monotone triangular transport maps. Foundations of Computational Mathemalics, 24(6):2063 2108, 2024a.
Baptista, R., Morrison, R.,, Zahm, O., and Marzouk, Y. Learning non- Gaussian graphical models via Hessian scores and triangular transport. Journal of Machine Learning Research, 25(85):1-46, 2024b. URL http://jmlr org/ pptprp//22//1110222222222 html.
Bartlett, P. L. and Maass, W. and neural networks, pp. 1188 1192, 2003.
dimension of neural nets. The handbook of brain theory
Bartlett, P. L. and Mendelson, S. Rademacher and Gaussian complexities: Risk bounds and structural results. Journal of Machine Learning Research, (:Noow)463--44 2002.
Bartlett, P. L., Bousquet, O., and Mendelson, S. Local Rademacher complexities. The Annals of Statistics, 33(4):1497-1153 2005.
Bartlett, P. T.,, Harvey, N., Liaw, C.,, and Mehrabian, A. Nearly-tight, VC-dimension and pseudodimension bounds for piecewise linear neural networks. Journal of Machine Learning Research, 20(63):1 17, 2019. Berti, P., Dreassi, E., and Rigo, P. Compatibility results for conditional distributions. Journal of Multivariate Analysis, 125:190 203, 2014.
Besag, J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society: Series B (Methodological), 36(2):192-22) 1974.
Billio, M., Getmansky, M., Lo, A. W., and Pelizzon, L. Econometric measures of connectedness and systemic risk in the finance and insurance sectors Journal of Financial Economics, 104(3):535-55) 2012.
Bühlmann, P., Peters, J., and Ernest, J. CAM: Causal additive models, high-dimensional order search and penalized regression. The Annals of Statistics, (2(6):252622535 2014.
Cai, T. T., Li, H., Liu, W., and Xie, J. Joint estimation of multiple high dishnsimensi precision matrices. Statistica Sinica, 26(2):445, 2016.
Cowell, R. G., Dawid, A. P., Lauritzen, S. L., and Spiegelhalter, D. J. Building and using probabilistic networks. Probabilistic Networks and Expert Systems, pp. 25-41, 1999.
Dasgupta, A., Murgoitio-Esannd; J., Ray, D., and Oberai, A. Conditional score-based generative models for solving physics- based inverse problems. In NewrIPS 2023 Workshop on Deep Learning and Inverse Problems, 2023. URL https: /tttps://ppenee
Diebold, F. X. and Yilmaz, K. On the network topology of variance decompositions: Measuring the con- nectedness of financial firms. Journal of Econometrics, 18991):119-134, 2014.
Dudley, R. M. VN Sudakov's work on expected suprema of Gaussian processes. In High Dimensional Probability VII: The Cargèse Volume, pp. 37-43. Springer, 2016.
Fan, J. and Zhang, W. Statistical methods with varying cocfficient models. Statistics and its Interface, 1 (1):179, 2008.
Friodman, J., Hastic, T., and Tibshirani, R. Sparse inverse covariance estimation with the graphical lasso. Bioslatistics, 9(3):432 441, 2008.
Greene, A. S., Gao, S., Scheinost, D., and Constable, R. Task-induced brain state manipulation improves prediction of individual traits. Nature Communications, 9(1):2807, 2018.
15

---

Published in Transactions on Machine Learning Research (05/2025)
Guo, J., Levina, E., Michailidis, G., and Zhu, J. Joint estimation of multiple graphical models. Biometrika, 98(1):1-15, 15, 2011.
Haussler, D. Sphere packing numbers for subsets of the boolean n-cube with bounded Vapnik-Chervonenkis dimension. Journal of Combinatorial Theory, Series A, 69(2):217-233) 1995.
Hyvärinen; A. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(),, 2005.
Khavari, B. and Rabusscau, G. Lower and upper bounds the pseudo-dimension of tensor network models. Advances in Newral Information Processing Systems, 34:10931 10943, 2021.
Kim, S., Cho, H., Lee, D., and Webster, M. J. Association between SNPs and gene expression in multiple regions of the human brain. Translational Psychialry, 2(5)=e1 el 13, 2012.
Klochkov, Y. and Zhivotovskiy, N. Stability and deviation optimal risk bounds with convergence rate o(//n). Advances in Neural Information Processing Systems, 34:5065-50777 2021.
Kohler, M. and Langer, S. On the rale of convergence of fully connected deep neural network regression estimates. The Annals of Statistics, 49(4)22231 2249, 2021.
Kolar, M., Parikh, A. P., and Xing, F. P. On sparse nonparametric conditional covariance selection. In Pro- ceedings of the 27th International Conference on International Conference on Machine Learning (ICML),
pp. 559-566, 2010.
Lafferty, J., McCallum, A., Pereira, F., et al. Conditional random fields: Probabilistic models for segmenting and labeling sequence data. In Proceedings of the 18th International Conference on Machine Learning (ICML), 2001.
Lauritzen, S. L. Graphical Models, volume 17. Clarendon Press, 1996.
Lee, K.-Y., Ji, D., Li, L., Constable, T., and Zhao, H. Conditional functional graphical models. Journal of the American Statistical Association, 118(541):257-271 2023.
Liu, H., Lafferty, J., and Wasserman, L. The nonparanormal: semiparametric estimation of high dimensional undirocted graphs. Journal of Machine Learning Research, 10(10), 2009.
Liu, H., Chen, X., Wasserman, L.., and Lafferty, J. Graph-valued regrossion. Advances in Neural Information Processing Systems, 23, 2010.
Liu, H., Han, F., Yuan, M., Lafferty, J., and Wasserman, L. High-dimensions semiparametric Gaussian copula graphical models. The Annals of Slatislics, 40(4):2293 2326, 2012. doi: 10.1214/12-ASS100
URL /ttpss://oo.. http:://dd..o/////// 10. 2214/12-A0S1037
Loh, P.-.. and Bühlmann, P. High-dimensional learning of linear causal networks via inverse covariance estimation. Journal of Machine Learning Research, 15(1):3065 3105, 2014.
Lounici, K., Pontil, M., Tsybakov, A. B., and van de Geer, S. Taking advantage of sparsity in multi-task learning. Proceedings of the 22nd Annual Conference Learning Theory (COLT), 2009.
Marcinkevičs, R. and Vogt, J. E. Interpretable models for Granger causality using self-explaining neural networks. In International Conference on Learning Representations (ICLR), 2021.
Meinshausen, N. and Bühlmann, P. High-dimensional graphs and variable selection with the Lasso. The Annals of Statistics, pp. 1436-1462, 2006.
Murphy, K. and Fox, D. Towards consensus regarding global signal regression for resting state functional connectivity MRI. Neuroimage, 154:169 173, 2017.
Ni, Y., Stingo, F. C., and Baladandayyuthapani, V. Bayesian covariate-dependent gaussian graphical models with varying structure. Journal of Machine Learning Research, 23(242):1 29, 2022.
Niu, Y., Ni, Y., Pati, D., and Mallick, B. K. aovsriateeassisteed bayesian graph learning for heterogeneous data. Journal of the American Statistical Association, 15977)1)11985 1999, 2024.
Ohtani, T., Nestor, P. G., Bouix, S., Saito, Y., Hosokawa, T., and Kubicki, M. Medial frontal white and gray matter contributions to general intelligence. PLoS One, (112): 2014.
Olver, F. W., Lozier, D. W., Boisvert, R. F., and Clark, C. W. NIST Handbook of Mathematical Panctions. Cambridge University Press, 2010.
16

---

Published in Transactions on Machine Learning Research (05/2025)
Pollard, D. Empirical Processes, volume 2. Institute of Mathematical Statistics, 1990.
Qiao, X., Guo, S., and James, G. M. Functional graphical models. Journal of the American Statistical Association, (14(525):211 222, 2019.
Ravikumar, P., Lafferty, J., Liu, H., and Wasserman, L. Sparse additive models. Journal of the Royal Statistienl Society Series B: Statistioal Methodology, (155:1009-10333 2009.
Song, M., Zhou, Y., Li, J., Liu, Y., Tian, L., Yu, C., and Jiang, Brain spontaneous functional connectivity and intelligence. Neuroimage,
2008.
Thompson, R., Bonilla, F. V., and Kohn, R. Contextual directed acyelic graphs. In International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 2872 2880. PMLR, 2024.
Vapnik, V. and Chervonenkis, A. Y. On the uniform convergence of relative frequencies of events to their probabilities. Theory of Probability & Its Applications, 16(2):264 280, 1971.
Wainwright, M. J., Jordan, M. I., et al. Graphical models, exponential familics, and variational inference. Foundalions and Trendso in Machine Learning, (1-2):1 11-2)):3305, 2008.
Wainwright, M. Sharp thresholds for noisy and high-dimensional recovery of sparsity using fl-constrained quadratic programming. TEEE Transactions Information Theory, 55(5):2183 0222183-22000 2202, 2009.
Yang, E., Baker, Y., Ravikumar, P., Allen, G., and Liu, Z. Mixed graphical models via exponential families. In Proceedings of the 17th International Conference on Artificial Intelligence and Statistics (AISTATS), volume 33, pp. 1042 1050. PMLR, 2014.
Yang, E., Ravikumar, P., Allen, G. I., and Liu, Z. Graphical models via univariate exponential family distributions. Journal of Machine Learning Research, 16(1):3813-335 3847,, 2015.
Zeng, Z., Li, M., and Vannucci, M. Bayesian covariate-dependent graph learning with dual group spike- and-slab prior. arXiv preprint arXio:2409.1144 2024.
Zhang, J. and Li, Y. High-dimensional Gaussian graphical regression models with covariates. Journal of the American Statistical Association, 118((53):2088 2100, 2023.
Zhao, T., Liu, H., Roeder, K., Lafferty, J., and Wasserman, The huge package for tigh-dimensional undirected graph estimation in R. Journal of Machine Learning Research, 13(1):1099-1022, 2012.
Zheng, Y., Ng, I., Fan, Y., and Zhang, K. Generalized precision matrix for scalable estimation of nonpara- metric Markov networks. In The Eleventh International Conference on Learning Representations (ICLR), 2023. URL https: //penreviieee http:///rrriiieeeeetourrrr///aaa
Zhou, F., He, K., and Ni, Y. Causal discovery with heterogeneous observational data. In Uncertainty in Artificial Intelligence, pp. 2383 2393. PMT.R, 2022.
Zhou, S., Lafferty, J., and Wasserman, L. Time varying undirected graphs. Machine Learning, 80:295 319, 2010.
17

---

Published in Transactions on Machine Learning Research (05/2025)
A Technical Appendix
A.1 Technical Definitions and A Useful Lemma
For notational convenience, the risk function is abbreviated to = =- トー el [Ett Pyu (z)(2))<=5)) for the remainder of the presentation. Let H == {h : y + R be a family of mcasurable functions, and {= E {-1,+11] collection of i.i.d Rademacher random variables. The Rademacher Complexity and Rademacher Auerage (Bartlett & Mendelson, 2002) are defined as:
n Σ o h(y"))
n sup - LEEH Ti 1=1
R, H -
sup = bEA TI
RH E fo
Σ o'l())))
Further, we define the empirical average and population expectation of h, as
Fnh = 1 Σ h(y)
and Ph: :h:=E1:0 Eyh(y))
In particular, we view y ()) and h(y) := (l (y) C(3; =) throughout our analysis. We further define the Local Rademacher Complerity and Local Rademacher Average with radius r as Rr H, P,P,, = rf and Ri E H. P12 < sr. The star hull of set of functions FF is defined as
+F == fa 10,
The following concepts are used in the proofs of the main results: Definition 1 (Lz-Covering Number). Let S, == {8 Vi be a set of points with y' E V. A set u C R" is an 6-cover w.r.t Lz-norm of F on { if vB E F, 3u EU, s.l. V M B(y')) S E, where is the i-th coordinate of u. The covering number N, S.) with L2-norm of F on S, is
minf[4]: U is E-cover of F on Sn).
and the covering number of F with Lorm of size n is N2, F,, n) == sy N2(E,F, Definition 2 (VC-dimension (Vapnik & Chervonenkis, 1971)). The VC-dimension dvc(H) of a hypothesis class H = {h: y + [1,-1}} the largest cardinality of the set S C y such that for all subsets S, denoted by S, 3h H:
(1 3
{{- -1 f
f(v) =
Definition 3 (Pseudo-dimension (Pollard, 1990)). The Pseudo-dimension dp(H) of real-valued hypothesis class H (h:3) [a, b]} is the VC-dimension of the hypothesis class
H = (h: V (-1,1) My t) = sign(h(y) - 1), E H).
Definition 4 ((m, C)-smoothness (Kohler & Langer, 2021)). Let m: ++ for some E No and A function f(+) [R9 + R is called (m, -sssmooth if for every (a1, a., ag) with M ag t, the partial derivative 0f 1/(00=222))) 21, .... 20)2g)( exists and satisfies:
0. (z) (001 21, goazg)
05
(g)
S C -
(0 goa Zg)
The following Lemma. is established which will be used in the proof of Theorem Lemma 1 (Verification of the Bernstein Condition). Under Assumptions 1 and the following inequality holds:
L L' F ((3) C(3)
F (c)) [(3))
2x
18

---

Published in Transactions on Machine Learning Research (05/2025)
Proof.
E Ee,=(((())
((B)) $ La, L2
Σ Kej
B3u (2) - Emees Pir (2)(() (x)
1=1
1.2 p
Σ
En,= |e
(EEnss) (2)( -e( {(Ensh Ti (x) )
Byk
2a
L
((3)
((3)- L(3)
20
A.2 Verification of Assumptions 1 and For the Mean Squared Error Loss
Note that the mean squared error (MSE) loss function (ie., E(a, b) (a-b)²) is used to recover the structure of the graphical model (sce (3)). Next, we establish that it satisfies Assumptions and 2.
Recall that we assumed that the sy's are uniformly bounded; e-g. without loss of generality we consider || I 1. Further, let (-x(-) E and F belongs to family of uniformly bounded functions: [|B())]l2 S 1.
Next, we establish that the MSE loss 1-strongly convex and 4-Lipschitz. Note that being quadratic function is trivially 1-strongly convex. The Lipschitz continuity claim follows from:
(G1 x) - ((f2, (31 (2) + 2x) S 4 T2/-
Next, we show that Assumption 2 holds; namely, Ex,> Ex,=1E)=1 (((1,, x x-1) 23 831 = 0. Recall that the posited DGP is:
8 2j
(z) + Ezz
Then, the following calculation shows that Assumption 2 holds:
Σ (13 a))) (87, x xr- o5))
( x. x.))
= F,,zz
j=l
F F
Σ (B, B7. x x-))
1=1
=0 if 18-88
A.3 Proof of Theorem 1
For ease and conciseness of the presentation, the hypothesis class H and its members defined the excess risk, the main quantity of interest in this analysis.
Proof. We start by defining the following function class:
= Ao (C L(3) C(B") : 3 F
Correspondingly, we denote h(x, x) == Ac, AC(3: xx,), which a composite function B from Fpx(9-1); loss function L and the difference operation A. It could be casily verified that y h E H. h(a, 2) e [0. In addition, h(-) also satisfies the Bernstein Condition by Lemma 1. On the other hand, directly bounding the Pseudo-dimension of H is trivial. To that end, we leverage the Lipschitz property of
19

---

Published in Transactions on Machine Learning Research (05/2025)
and the covering number of F to construct an upper-bound the covering number of H. We will focus analyzing the Local Rademacher Average ER, {h E H p
S r):
t E H. pnh Sr) = E
sup NER.Pa
1 Σ
o''( 25)
sup BEF P, ~S 4.9
! Σ 0
1=1
i=1
Next, we bound PA c.B PnA 2. To invoke Theorem 3.3 in Bartlett et al. (2005), we need to find subroot function r(r) such that
c.8
2L LR, th E H: [[h"] r)
7()) 2
a
By Lemma 3.4 from Bartlett et al. (2005), it suffices to choose:
*(r*) 207,2
ESt, (*H, sr) 11 log
a
n
with fixed point y" r(r") denoted r*.
The following analysis largely follows from the proof in Corollary 3.7 in Bartlett et al. (2005). Since A c.3 is uniformly bounded by 1, for any r 2 Tr), Corollary 2.2 in Bartlett et al. (2005) implies that with probability at least 1, {h *H. 1hh2 I < S (h E *H P2 < 2r}. Let & == { E #H PA S r {h *H
P, h,h? < 2r}, then the following holds:
ER,, (HH, Ph2 1r) S P[ HH Ph S PIE PH? < <rH8
s E[ h2
2r + n2
Since r* is the fixed point of a sub-root function, namely rr = T(r"), by Lemma 3.2 in Bartlett et al. (2005), r" satisfies the following
20L r" & ED; (HH, & 2r") 11 log n 20
(10)
n
where the Lipschitz constant L and the strong convexity parameter show due to the Bernstein condition.
Next, we leverage Dudley's chaining bound (Dudley, 2016) to upper bound (HH, < 2r"}, using the integral of covering number. Specifically, by applying the chaining bound, it follows from Theorom B.7 (Bartlett al., 2005) that
Es, PR (R,,,,,,,,,, (HH, < 2r*)]
const
Vn Es.. /
VlogN2 *H,
(11)
where const represents some universal constant. Next, we bound the covering number N2(e, S,,S.) by (ele/LL,FFx((p11 S.). We show that for all Sn, any -cover of Fpx(e-1) is é-cover of H, which implies that N2(e, S,) S N2(e/ Lp, Fp Specifically, let Usz C [0, 1|" be an e-cover of F on S, so that for all 3it E F, 3 3 Ha E Uik so that
Σ (3zx(z)) - w//))² SE. dE|n|
n
Further, let U C Raxpx(p-1), B(z) 2 RPx(1-1) where B(z) E Fpx(p-1) is a family of p (p 1) joint functions, with each element Bi E F.j,k E [p].] & k. We say U is an e-cover of Fpx(g-1) {z 1= if v/ E Fpx(p-1); Eu E U s.t.
Σ
Σ
(ug, Br(()))) Se.
np(p
ie [m) je(n) Agj
20

---

Published in Transactions on Machine Learning Research (05/2025)
Clearly, let Uzk be any collection of arbitrary p x (p- - 1) E-covers of F, the Cartesian product of J E 1 + forms an E-cover of (x-11-)); which implics that | S [M pp Thus we have (e, Px(e- Sp) S (E, F, S.))" Next, we show that given any ip -cover of F (p-1) on Sn, denoted as 4t one can construct =1,))))))))))))
V = (v = (w',-. ... vny E Rm|0 = D (e((u),, x² -))))
[n], u U3, which is an
E-cover for H. i.., for all h 3.3 V that v 1/2 EiE 22) ud) S
2 e(18) mb) A cas 21)
((uu) x. m.)
Σ
sE[n]
P
x
sElPl
2
((lu), x² =)))))))))))
e(133, x² -))(zz))
P
P JE(p]
sE[n]
jE(p]
[uu G.. x²
n2p2
kc(m)
E
L2
((u)) - 9g. (ل.
np
SEIn| jElp
(3)
ELp.
In the above derivation, (1) leverages the fact that.
Ac,g( 2²)
Σ e (8)) a- .²)))
-1 Σ
(183
p selpi
jEiPl
(2) is by Lipschitz continuity, and (3) is by the definition of the covering number and the fact that |z]]| S
1, vk [p]. Combining the above incquality with Corollary 3.7 from Bartlett et al. (2005), we have
(121+1)
S,, log
N₂
2 H, Sm
< log { { N/ 8LP Fp Sm
(
(1²
s p² log No F,
(12)+1)}
(11+1
8LP
Next, we bound const E 22 Vlgg vn number (Haussler, 1995) we have:
VlogNN(E, N2(,, *H, dne from (11). Note that by Haussler's bound on the covering
log
s const - log ( F, Fr
E
Lp Sn,
8T.'
and therefore
const
/2r-
const
[o V2r*
Nogg (
(12)+
F.S.
V Nag,, *H, Sn)de=
Ess.
2 H, Sn
Vn
Vn
const
P
Es
//2
log.V2
Fr ) ( [ + 1
) de
Vn
8T.p"
dr(F) pva
Tpp log
dr(F)r* log(T./r*)
const PV
< const PV
n
E
n
(F)
dp(F)r*
const PV
n² +
(12)
21

---

Published in Transactions on Machine Learning Research (05/2025)
where in (12) refers to Euler's constant and const represents some universal constant that may change from line to line in the above derivation. The inequality in (12) comes from the fact that r loogr is monotone increasing function for r* e: in the case where r* S ad n. loge(n/ed] S we have dr(F)r' log(1/r') S di; in
m
the case where r* 2 ed (n/ed// we have dp(F)r" log(1/r") X dp(F)r" log(n/edp(F)). Combining these two
cases yields inequality (12).
Together with (10) one can solve for
dd (og(F))og((a uLp
(p())
a²n
By Theorem 3.3 in Bartlett et al. (2005), have that for all B px with probability at least 1 6:
R(3) - R(8)) S RR RR. + dp(F)
p2pp² log(nLp) log
an
Rw Rx) - (8*) & R(B) R(8") + dp(F), 222²² log(nLp) log
(13)
an
which gives
R(3) - R(3) < R3) - Rr (8)) + dr(F), L2pp² log(nLp) log
(14)
eXcesss crrOr
empirical excess error
generalization
Due to the fact that R(3) is the empirical risk minimizer, we have:
Rr R.) Rr (3*)
(1,(2') xe )) Σ
(1) 2 x0 x))
empirical excess error
(E(m)
j=1
5=1
(B8100) (2') xx. =)
Σ
8 (18](²)) 20 و
(15)
np
vE(x)
/=1
1=1
L
|| 1opt 37
(16)
np
(r(
j=1
<
L Expprox
approximatiio error
Equation (15) is due the fact that is the empirical risk minimizer within F Fpx(p-1) and thus R, RR, (gopt. Finally, Equation (16) is derived by the Lipshitzness of E(,,)).
A.4 Proof of Corollary 1
Proof. Let e : p²("d an
(Tp()) log (?), by (6) we have R(Ø) R((8") $ T. Eapprox(F) + e, which could
be expressed as:
E
Σ #=1
e (0)()) x Σ 1=1
(3)22),=-))0000 (197())) $ L. Expprox (F) + E.
By the a-strongly convexily assumption (Assumption 1), namely,
((a1bb) - E(a2; b) 2 E'(a; b) (a1 az) + a (a1 - az)²
2
22

---

Published in Transactions on Machine Learning Research (05/2025)
by setting 01 (8,(2)) x))) 02 === (z) xz)), ;)), b= z, and the optimality of Øj. we have:
5 + L. E-pprrxpro 2 Ex
Σ
(1)
- -E
l (83(2); x xc )]
j=1
E F,,
(((1,, x -)) 3z , = (B)(2) (z(2),, xz
=0 by optimality of 9
P
a
Σ
2
j=1
which implies that
(3)))) En 5== T | (8)(=) 83(3))1
E Σ (B)(=)
= EE
j=1
28
2L Eapprox
A.5 Proof of Theorem 2
The result aims to balance the function approximation crror (O(Eapprox) and the genoralization error (0(de) de)) )).
Proof. The major tool that we leverage is from Kohler & Langer (2021). particular, we invoke Theorem from Kohler & Langer (2021) with 1, p m, M e- 1/1m L ~ { 4/2m r (2e)²( (2 (m+t)q², we have that Vie [p].j E 3 gop F, 11 3ixllso S & thus Eapprox (F) s &. On the other hand, by Theorem 6 from Bartlett et al. (2019), we can bound the VC-dimension in terms of number of layers and neurons. In the case of deep neural network (Theorem 2b from Kohler & Langer (2021)):
jk
dp S const (²)*log(Hr)
(BeM)) m (m+t9) 4 1 m x²) log(Hr)
< const.
S const
(22 (log(22 M) + log(m) + log(g)).
It suffices to pick M = const
8- 1/2m and one can get that
dp const -g- 5- m2 log(E).
Plugging the above into (6), one has:
R( - R(8)) S + L2g-a/ m/ p² (oo((1/5 (og(nL)))))))))l)))))))))))) log (
(17)
an
To minimize the right-hand side of the above inequality, it suffices to choose { Lmmgppp log" (nLp) log(p) log(
(18)
am
Finally, Equation (8) could derived using a similar strong convexity argument as in Corollary to bridge the gap between excess risk and the quality of B: [=]])((2) (=))))a(()))
23

---

Published in Transactions on Machine Learning Research (05/2025)
A.6 Extension of Theorem 2 with An Additional Assumption on the Loss Function
In this section we provide an exlension of Theorem 2 under the following additional assumption:
Assumption 6. Assume the following holds in addition for the loss function E(:--), namely its gradient is also Lipschitz:
b(ar;b) - b(a2;b))== & (aaaba)a=a (a1 -2) + 1 - -a2)²
This is equivalent to assuming that the second derivalive of the loss function ("()) is upper bounded. Next, we state the result formally.
Theorem 3. Let F correspond to family of fully connected neural networks with ReLU activation func- tions, with number of layers H = g-9/2m and number of neurons = (2e) (m++))²² Then, under the assumptions of Theorem 1, logether with Assumption 4 with m S q and Assumption 6, by setting [m [mm'dpp d2 log (nLp) log(p) log(1/0) log(/a))
2m-4
(19)
Aan
the following holds with probability at least 1 -:
8:
1 P
x Xe
R(B) - ((B) S Ag" F-
(2) -
(20)
Proof. Using a similar argument as in the proof of Theorem 1, starting from (15):
RT) - (g") =
(8)(z)), x" m
Σ
(97 x²) x x
mp
kE|n]
j=l
3=1
(0000 (z'), x 5)
up sE/n|
j=1
RM. Rr (B*).
With incquality (13), the following holds with probability at least 1 - 8:
(BOPt) - (3') ( R(BoU)) RR(B) + dF(F) r,2p² log(nIp) log
on
Together with (14), we have:
R(3) - R((3)) S R( - R(8*++++++++ dr(F) L2pp log(nLp) log i
(21)
an
For the term R(gort) 1 R((3²) on the RHS, it satisfies
e (1890(( (x), 1)) - Σ 8 (8)( x
R(Bort) - R(B")- =
Eg,x x
npp
1
I
(1B 31 x- -)) 11) -8 (0)p0(z) (3) x
np
=0, by the optimality of 87
E
Σ
1(8 (2) - -)(2)) x;
2 np
19
1.80pt
171120
2np
[m(m]
1=1
<
Eapprox (F)2
2
24

---

Published in Transactions on Machine Learning Research (05/2025)
Plugging the above inequality into (21) we have:
()8)-R(03)) x A (F) + dr(F) L2pp log(nLp) log (7)
on
Let Eapprox (F) = &, by plugging inequality (17) into above inequality, we have:
R(0) - R(3") 2 5² + 5 5///" L²pp² log (Lp()L log
1 8
on
To minimize the RHS of above inequality, if suffices to choose
~ Lm log" (nLp) log(p) (1)) log(1/a) Xan
Remark 10. The choice of É in (19) implics that the R(3) -((()) converges at a rate of the order n 2m amm/2m++))) which matches the rate from Kohler & Langer (2021).
A.7 Proof of Corollary 2
From inequality (7) we have:
P
Σ
A,(=) (3)3()
12dp(F) 1² logg(IIp) log
L E Epproxx (F)
2 -2(B-8)21 -
(5)
+
al3 - 332
Using Assumption 5, we have (0)(2) (2(2), A,(z)(3 (3()))
1113(22) (2)2))11 for each
1, p and therefore:
1 (8)(()) S L2d p² log(n log 3 a²>(B 3)2
+ L Eapprox (F)
E.
a7 -322
On the other hand, by the fact that the margin between strong and weak edges is bounded away from zero, namely, (2 (0- B) 2 & > O, we can select a threshold r = n1 + (1 - n)ß and measure the edge-wise guarantee using the binary risk: 1=)3 2 +} * [1)) 2 $}.
Next, we establish the following incquality:
9 < -(xx)))))
min(??, (1 (1- -
(22)
The starting point is the following decomposition:
( = 11 () 2) (g)))) - -B)())))))))))
Term
21] (z) 2 8) (B)((() -8,(2))
Term II
To show inequality (22), it suffices to bound Term In case 1113 2 # 11135((2) 2 3}, we have either [))(2)) 2 7, 18185x < 3, or [3)x(=)) < 7, (=) 2 B. In both cases, we have min(n? -
(1- )) 3) <
(z)(2) - (=)), and hence
min(112 - () Bipe-BBp(n
B8 (1(11188(()) 2)
#1) + 2
) s ( - -3)(2))²
25

---

Published in Transactions on Machine Learning Research (05/2025)
The last expression implies that
p
min()²,(1-7)))6
ΣΣ kØ
113x(()) + +11))())1
51 - Σ> Σ
- -36(2))
y=1 Kgj
Consequently, it follows that
P 2 #1 2 ()
F
ΣΣ k=
[13x(())1
(F)
L. E.Eaproxx (F) 07 min (1p, (1 (B B)2"
min (iin ( 33)mm log(nLp) log {
B
Additional Details and Results for Synthetic Data Experiments
Additional illustration for the data generating process. Figure provides visualization for candidate skeletons y'' that are used in settings G1, G2, N1 and N2, and selected O's that are obtained as the convex combination of the candidates, depending the value of their corresponding zss.
doeletoo
sieltton M.
clazer H-1,
sample
(8) skeletons for 11 22, 2s for settings G1 and N1
(b) examples for 9''s whose skeleton depends 2
aoelatoo My
chater idel, sample Mde
clatter Meni, sample Idnd
chaster idv2, samcle
skeletons for 1. 2. Va, for settings G2 and N2
examples for 0¹ whose skeleton depends 2
Figure Pictorial illustration for settings G1, Left: candidate skeletons V1, V2, Vs. Right: different On's that obtained from candidate skeletons, with the exact mixing depending the values of the corresponding 2"s; their diagonals are suppressod for visualization purpose.
Figure 3a provides visualizations for the candidate skeletons B1, Bg whose corresponding graphs have tree structure and their convex combination (no longer corresponds to a tree); they serve the skeletons for the DAGs. Figure 3b shows the resulting G's after moralization. The color (red/blue) and shade respectively reflect the sign (positive/negative) /ngaatiee and the magnitude of the entries in the moralized graphs for the linear case where the exact values can be calculated.
Of particular note, in the special case where the DAG possesses tree structure, the moralized graph is equivalent to removing the direction of the edges in the DAG, and therefore the skeleton matrix can be obtained by adding the "transposed" entries to that of the DAG; e.g., the first two plots in Figure 3b. Otherwise, in addition to the "transposed" entries, there are additional "married" entries as result of connecting nodes having common children. As it can be seen from the rightmost plot, these "married" entries (the very faint ones in the middle of the heatmap) tend to be extremely weak in magnitude, and pose challenges for estimating them due the very low signal-to-noise. Finally, for the non-linear case, one can
Note that by definition of the DAG, nodes possess topological ordering and thus the corresponding matrix can be written as triaguuarr after reordering; here we showing the skeleton matrix after such reordering.
26

---

Published in Transactions on Machine Learning Research (05/2025)
custer del, somole Me.
clette (de2. sample kied
Guner de3 sample
(a) skeletons for B,, Bg and an example of their evx combination
(b) moralized graphs corresponding to the DAGs in 3a
Figure 3: Pictorial illustration for sctting D1. Left: candidate skeletons B1, B2 and example of their convex combination.
Right: Moralized graphs for the respective DAGs, with their diagonals suppressed for visualization purpose.
oblain values of the moralized graph approcimalelyy by considering linear approximation to the ʃi&" 's (recall that z) = Enepat) (z+) + te),, and they exhibit qualitatively similar patterns to the linear case.
Additional simulation results. Table 4 presents additional metrics for the DNN-based covariate- dependent graphical model estimation method. In particular, given the sparse nature of the underlying true graphs, we report. the F1 score (2 precision-recall and balanced accuracy (Snnatiiity+seecicctt) when the 2
precision-recall estimated graphs are thresholded at different levels. The result shows that for practical purposes, practi- tioners can effectively recover sparse skeleton that is close to the truth by applying reasonably small thresholding.
Table 4: score (F1) and Balance Accuracy (BA) for the estimated graphs by DNN-GCM at different thresholding levels.
G1
G2
N1
N2
D1
D2
threshold
F1
BA
BA
F1
BA
F1
0.010
0.23
0.71
0.15
0.84
0.23
0.71
0.15
0.84
0,23
0.76
0.24
0.76
0.025
0.49
0.91
0.72
0.98
0.50
0.91
0.71
0.96
0.56
0.91
0.63
0.85
0.050
0.83
0.98
0.70
0.81
0.89
0.99
0.65
0.77
0.75
0.89
0.56
0.73
0.075
0.91
0.98
0.57
0.71
0.95
0.98
0.56
0.69
0.70
0.84
0.32
0.60
0.100
0.94
0.97
0.48
0.66
0.93
0.96
0.41
0.63
0.65
0.79
0.02
0.51
Note: reported metrics are first averaged across test samples for single experiment, then averaged across experiments 5 data replicates.
Table 5 presents additional results for settings D1 and D2, where instead of evaluating the performance against the true moralized graph O³, we evaluate against a "pseudo" moralized graph Ô' where we treat odges that are present due to married nodes as they were non-existent. As previously mentioned, these edges typically exhibit very small magnitude and hence are difficult to recover. We expect the performance improve slightly if the comparison were done against the Gi's, which is indeed the case: AUROC improves by around 3% and AUPRC improves by around 5-7% (see, e.g., the last two rows in Table 3 for comparison).
Table 5: Performance evaluation for DAG-based settings against the pseudo moralized graph without edges from married nodes.
DNN CGM
RcgGMM
glasso est. by cluster
nodewise Lasso- est. cluster
AUROC
AUPRC
AUROC
AUPRC
AUROC
AUPRC
AUROC
AUPRC 0.89 (0.020)
D11
0.98 (0.003)
0.81 (0.028)
0.97 (0.0U6)
0.81 (0.028)
0.97 (0.008)
0.50 (0.007)
0.99 (0.002)
D2
0.96 (0.013)
0.74 (0.032)
0.91 (0.015)
0.65 (0.023)
0.94 (0.009)
0.44 (0.002)
0.92 (0.002)
0.52 (0.110)
Remark 11. For DNN-based method, the performance improve when training is done with larger sample size. In particular, for challenging settings such as D1 and D2, for estimates obtained from model trained using 30,000 samples, when they are evaluated against the true graph O¹, AUROC is given by 0.97 (0.004) and 0.95 (0.007) respectively for D1 and D2, and AUPRC given by 0.91 (0.004) and 0.75 (0.015). When the estimates are evaluated against the pscudo moral graph O' (so that extremely weak entries from married nodes are not counted toward the skeleton), AUROC reaches around 0.99 for both settings and AUPRC reaches 0.95 and 0.82, resp.
27

---

Published in Transactions on Machine Learning Research (05/2025)
Finally, Table 6 shows the metrics for glasso and nodewise Lasso when they perform estimation on the full set of samples; note that in the absence of covariate- estimation methods, this would be how the two methods perform in practice since the partitioning would not be known apriori.
Table 6: Evaluation for glasso and nodewise Lasso on full samples without partitioned by clusters.
G2
N2
D1
D2
D1 pseudo moralized graph
D2 peeudo moralized graph
AUROC
0.95
0.96
0.91
0.89
0.95
0.94
glasss
AUPRC
0.17
0.46
0.47
0.45
0.39
0.37
0.43
0.11
nodewire Lasso
AUROC
0.97
0.98
0.97
0.97
0.92
0.88
0.97
0.94
AUPRC
0.66
0.67
0.66
0.47
0.55
0.43
0.62
0.53
Note: values correspoud to the average over experiments on 5 data replicates. The last 2 columus show the evaluation against paeudo moralized graphs under DAG settings.
Compared with the results in Table 3, nodewise Lasso exhibits material improvement when the dependency the covariate is linear/milan non-linear, whereas glasso seems more susceptible to the mis-specification, induced by the varying magnitude across samples, as manifested by a much smaller improvement in perfor- mance in the case where estimation is conducted on partitioned samples.
C Additional Real Data Experiments
In this section, we present results obtained from applying the DNN-based covariale-dependent graphical model to finance dataset involving S&P 100 stocks. Such model estimates the partial correlation across stocks, while conditioning on covariates that correspond to the broad market condition.
The S&P 100 constituent dataset. We examine the coeeeceonnets of the S&P 100 Index con- stituent stocks under different market conditions. These constituent stocks correspond to 100 major blue chip companics in the United States and span various industry sectors. In the sequel, stocks and tickers may be used interchangeably, and they both correspond to the nodes of the network of interest.
Data collection and preprocessing. We collect daily stock return (calculated based on adjusted close price) data for those that are components of the S&P 100 Index as of 2023-12-29. At the initial data gathering stage, the first historical date is set to 2000-01-02, and the last to 2023-12-29. We require all tickers to have valid data on the first historical date, and the set of tickers are further filtered ensure this; tickers such as ABBV, GM, GOOG, META etc are therefore excluded from subsequent analyses. After this filtering step, the remaining set encompasses 79 tickers and thus the size of the network p 79.
To obtain samples for x, we consider beta-adjusted SPX residual returns, so that the market information is profiled out from the mean structure. Concretely, as a data- preprocessing step, let m² be the S&PP 500 Index return on day i, which is a proxy for the market return of US large-cap companies; for each ticker indexed by j, let 85 be its raw return. We first estimate its beta 33 by regressing its return against that of the S&P 500 Index, over a lookback window of 252 days; the beta-adjusted residual return is then given by 2² == # - -3mm.. This gives rise to {aff. 1,..j=1..... 1, .P) residual return data (namely, the 's's) effectively starts from 2001-01-02 and there are total number of 5785 observations. Note that the processed data exhibit de minimis serial correlation and the æ can be regarded i.i.d. for practical purposes. For 2, we simply consider set of variables that can reflect the high-level market condition of the day, and it encompasses the returns of the S&P 500 Index, the Nasdaq Composite Index, respectively, and VIX. As such, q 3. Data is further split into train/val/test periods that respectively span 2001-2017, 2018-2019, 2020 onwards. Note that such split is more for the purpose of demonstrating that the trained model is capable of producing interpretable results when applied to unseen data (c.g,, COVID period); however, results are presented for all periods to demonstrate the patterns the proposed method uncover.
for all tickers. After this pre-processing step, the
Results. Figure shows various graph connectivity metrics over time (left) and estimated graphs selected dates (right), after thresholding. The observation is 3-fold: (1) the estimated networks exhibit
28

---

Published in Transactions on Machine Learning Research (05/2025)
significantly higher connectivity during volatile market conditions". (2) Under normal market conditions, intra-sector connectivity level is higher, as manifested by the corresponding heatmap showing a block- diagonal pattern; however, in adverse conditions, such a pattern is diluted due to increased connectivity across the board. (3) Overall, the graph dynamics are slow- movvnng in the sense that the estimated graphs usually show high concordance for dates that are not far apart.
66-
(a) graph connectivity motrics (top) and (b) heatmaps for estimated graphs on selected high-vix (left, 2020-03-23) Jow-vix level (bottom)
(right, 2023-03-23) days.
Figure 4: Left panel: varions connectivity metrics of the estimated graphs (top) and the corresponding VIX level (hottom), from the beginning of 2001 to the end of 2023. Notably events that significantly increased market volatility have been marked. Right panel: heatmaps of the estimated graphs on representative dates that respectively have high and low-VIX, with red cells indicating positive partial correlation and blue cells negative ones.
Figure 5 presents the average partial correlation graph of high and normal VIX days in the test period, respectively, where high VIX days encompass dates on which VIX exccoded the 90% percentile of the period in question and normal VIX ones correspond to those that were below median.
83
(a) Averaged graph for high- VIX days
(b) Averaged graph for normal- VIX days
Figure Averaged partial correlation graphs for high-VIX (left) and normal-VI (right) days, with red cells indicating positive
partial correlation and blue colls negative oncs.
After the averaging step, the increased connectivity during the high-VIX period is not pronounced that on the extremes (e.g, 2020-03-23) where the intra-sector pattern no longer stands out. Nonetheless, compared with the average graph from the normal- Mamal--IX&a days, the connectivity is still at an elevated level, as manifested by the additional inter-sector connections while the intra-sector ones stay comparable.
Similar patterns have been reported in existing financial cconometrics literature (c.g., Billio al., 2012; Diebold & Yılmaz, 2014), albeit models/methods methods therein examine different quantities rather than the precision matrix. Note also that these methods do not support sample-specific estimates and thus their results were obtained through rolling window analyses.
29

---

Published in Transactions on Machine Learning Research (05/2025)
Finally, note that refrain from dictating any case of the estimated graphs, which highly depends the specific context and application.
Graphical Model Preliminaries
D.1 Decomposition of Gaussian Multivariate Distribution into Clique-wise Potential Functions According to Graph G.
Consider a Gaussian random vector æ E RP with mean vector μ == 0 and covariance matrix E assumed to be positive definite. Its joint distribution function is given by
1 E 'z). P(x) (dettc)) exp(
Then,
log(2)) C +
Σ j.k=1
and hence the potential functions over pairwise cliques correspond to
=
they are functions of E which is the parameter of the specific Gaussian distribution under consideration. Further, it can be easily seen that zj TL 23 (-.A,)): if and only if Σ5 = 0. Hence, to estimate Gaus- sian graphical model from data, it suffices to estimate the inverse covariance matrix and all the pairwise conditional independent relationships encoded by a graph correspond to the zero elements of Σー
It has be further shown that the Gaussian graphical model can also be estimated based on regression tech- niques; namely, by estimating the regression coefficients of the following model (Meinshausen & Bühlmann, 2006), procedure referred to neighborhood selection:
Jj
Byrtn +
k= k=1,
where 35k := : /
D.2 An Exponential Family based node-conditional Graphical Model
Consider a random vector æ E RR whose joint distribution is P(x)). Further, consider the following collection of node conditional distributions from the univoriate exponential family, where the neighbors of node j according to an underlying graph G denoted by No(l)):
p(zzlo-;) x exp (0xx))(0+ +
Σ Oxng(xx) +
Σ
(jee((a)g(z))
(23)
kENG()
k,EENoU)
Σ
1 Ojmy me II g(xx)))
where the multivariate canonical parameter 0 is defined as linear combinations of up to C-order products of positive univariate functions g(xx) of neighboring nodes to j according to the graph G.
An application of the Hammersoley-Cliffor theorem (see Scction D.3) and some algobra (Yang al., 2015) show that the collection of the node conditional distributions in (23) define proper graphical model with respect to G.
30

---

Published in Transactions on Machine Learning Research (05/2025)
D.3 Markov Properties of Graphical Models
The conditional independence/ dedpreddeeeeeeeeeeeneeeeeeee relationships defined by graphical model are a consequence of the Markov properties encoded the graph G. Specifically, the probability distribution satisfics the pairwise Markov properly with respect to G, if æy JL I& ==(jA) for all [k,j] & E i.e. variables æj and we are conditional independent given all the other variables, if they set E of G. It satisfies the local Markov property with respect to G, if Iy TL (()(()c ( i.c., the conditional distribution of variable Ij given all its neighbors is independent of any other nodes in G. Finally, satisfies the global Markov properly with respect to G, if for any subsets A, B, C of V such that C separates A and B (i.e., every path between a node in A and a node in B contains a node in C`), the following relationship holds: EA IL 2B sc. It can be shown that the global Markov property implics the local one and both imply the pairwise one.
not connected according to the edge
Next, we elaborate the relationship between the joint probability distribution F(x) of the random vector æ and the underlying graph G regarding conditional crependence/depenndeee between the graph G and the probability distribution PP(x) comes through the concept of graph factorization. Specifically, let C(G) denote the set of all cliques of G. Then, P(x) factorizes with respect to G, if it can written as ((x) = Ž IICEC(a) Φc(#c), with dc > for all C E C(G); i.e,, the joint distribution of w can be written as a product of positive functions (called potential functions in the literature) that depend only on a subset of the random variables in the clique. However, the concept of graph factorization does seappeendenne/eepeeeeeeeee
relationships. The connection
relationships between subsets of the random veclor x. Such relationships are associated through the various Markov properties of G, as outlined above (see also Lauritzen, 1996, for indepth prosentation). Hence, if P(xx) factorizes over G, and G possesses Markov property, then conditional independence relationships for subsets of æ present under P. The reverse relationship, namely that if P satisfies a Markov property with respect to G, then it also factorizes with respect to G, is established by the famous Hammerseley-Clifforr theorem, provided that P(x) possesses a positive and continuous density (Besag, 1974).
not imply per se conditional
E Notes on Implementation
Implementation of the neural network. For all synthetic and real data experiments, we use fully connected MLPs (i.c., stacked modules consisting of Linear-ReLU-Droparr layers) with residual connections as the underlying architecture to parameterize the (-)(s)'s. In particular, we adopt the implementation scheme as mentioned in Remark 1, where we use a single neural network that takes z as the input (i.e. input dim is q); but with multi-dimenssonal output head that produces outputs for all Byx (,,()),,kk j. (i.e., output dim is p(p 1)).
[p], # j
MLP
hilden rep
Pia 81 8h
Pa.1
Pas ... 8ap-1 8:00
MLP
Ppu
Ba2 Pba
Bap-1
cogy
impat
COMENE
Outpus
size
slon
plp
Figure 6: Diagram for the A(=): RR RP(p-1) network; an MLP block consists of Linear-ReLLU--Dppppoe layers.
Overall, we observe that the above-mentioned implementation scheme is fairly robust to the choice of hyper- parameters; in addition, it shows material improvement in performance when compared with a single MLP block without the residual connection, when q becomes moderately large. Other "intermediate" sharing- backbone schemes have been experimented, yet the performance generally inferior and more sensitive During training, we fix batch size 512 with gradient clipped at 1.0, while
the choice of other hyperparameters are tuned over a hyper-grid, with the best set determined based the performance
31

---

Published in Transactions on Machine Learning Research (05/2025)
validation set of 1000 samples, with MSE as the criterion. All experiments done NVIDIA RTX A5000 GPU.
Table 7: hyper-parameters for the MLPs and model training
hiddon layer size /dropout
learning race
scheduler type schoduler stopsize
decay
epochs
G1
128.
64).
128
0.3
0.0005
StepLR
20/0.25
50
G2
128,
64
[28
0.3
0.0005
StepLR
20/0.25
80
N1
1288,
128
0.0005
StepLR
20/0.25
N2
[128, 64.
[128
0.3
0.0005
StepLR
20 20/0.25
80
[64,
32]
StepLit
20/0.25
D2
[64,
32}.
[64]
/ 0.1
0.0005
StepLR
20/0.25
80
In regards to the benchmarking methods, for glasso and nodewise Lasso, we rely the implementation in R package huge!s 18 (Zhao et al., 2012); for RegGMM, we did not find any publicly available packages/mode modules and thus implemented the method ourselves leveraging PyTorch, where (z)) effectively reduces to a network with a single linear layer.
Data and code availability. The resting-state fMRI scans dataset was provided by the authors of Lee et al. (2023). The S&P 100 Index constituents dataset can be collected from YahoolFinance with the list of tickers corresponding to the constituents available through Wikipedia The code repository containing all the implementation is available at htt.thaww.coo/eeoggecccaaaa.aa.a depeedent-graahieeeeee
iPnttpe://finaaceey/aee www.aa..ce.aa.ooooooooooo
hhttpp://em.uuikipddis...
100#Components
32


Reviews:
[{"review_id": "yzhxQFbFCl", "reviewer": "Reviewer_e16v", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strong aspects:\n- Theoretical derivations, assumptions, and results (PAC-style guarantees and finite-sample error bounds) are rigorous and seem technically sound.\n- Rigorously combine deep neural network-based functional approximation with covariate-dependent graphical models, providing a clear theoretical foundation.\n- Reasonable experimentation with appropriate baselines (linear and static graphical models), clear illustration of scenarios where the proposed approach excels (particularly nonlinear dependencies).\n- Does not require Gaussianity; accommodates general non-Gaussian distributions using node-wise regression and score matching methods, substantially broadening applicability.\n- Clear real-world examples demonstrating interpretability, such as cognitive-state-dependent brain networks and finance market volatility-dependent stock correlations.\n\n### Weaker elements (minor, but worth attention):\n- The process of translating continuous estimated weights to a symmetric undirected graph is not fully clear.\n- Not clear if regularization was used in the experiments and what are the parameter tuning strategies (e.g., choice of sparsity penalty or threshold).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The submission introduces a deep neural network-based methodology for estimating covariate-dependent graphical models, where the structure of a graphical model (representing conditional independencies) can vary continuously with external covariates. The paper provides finite-sample statistical guarantees (PAC-learning type bounds), ensuring theoretical rigor and convergence rates under stated assumptions. Empirically, it demonstrates superior performance over existing methods in both synthetic scenarios and real-world datasets (neuroscience and finance).", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Despite the overall clarity, there are a few areas that could be improved or elaborated:\n\n1. Symmetry of the Graph: The model (1) treats each node’s neighborhood separately, and the neural network outputs $\\beta_{jk}(z)$ potentially for each ordered pair. It’s implied that in an undirected graph, $\\beta_{jk}(z)$ and $\\beta_{kj}(z)$ should both be zero or non-zero when an edge is absent or present. Unless I missed it, the authors did not explicitly state how they enforce or utilize this symmetry. In practice, they might take the union of edges predicted in each direction or average the two estimates. A brief note on this (e.g., “we treat an edge $(j,k)$ as present if either $\\beta_{jk}$ or $\\beta_{kj}$ exceeds the threshold” or similar) would clarify the estimation of an undirected skeleton from directed regression outputs. This is a minor point, as it’s understood from context and common practice in neighborhood selection, but mentioning it would remove any ambiguity.\n\n2. Regularization and Tuning: I don't think there is an explicit mention about the form of the regularizer used, if any, in the main text. I suspect the authors likely used an $\\ell_1$ penalty on the coefficients to encourage many $\\beta_{jk}(z)$ to be exactly zero. An equation for the full objective (loss + penalty) would be helpful. Also, while the experiments indicate the authors thresholded the learned weights at a certain level, it’s not fully explained how that threshold was chosen. Did they use a validation set to pick an optimal threshold or lambda for the penalty? Or did they simply examine a range of thresholds (as in Table 3) to show stability? A sentence on the model selection procedure (e.g., using cross-validation for any hyperparameters, or using BIC to select a sparsity level) would strengthen the presentation.\n\n3. Summarize briefly differences and similarities with prior approaches (linear-varying, kernel-based, and score-matching methods) in a concise comparative table in the related-work section. This would facilitate easier comparison and highlight clearly the novelty and advantages of the proposed approach.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "hZ8EKNafjo", "reviewer": "Reviewer_Ppck", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- This paper provides an interesting theoretical result.\nIn particular, I have never seen the usage of Dudley's integral in this context.\n\n- The presentation is generally clear.\n\nWeaknesses:\n\n- Though the authors provide experimental results to support their theoretical guarantees, the model considered in the paper may be oversimplified.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the following problem.\nSuppose $G$ is a graph with $p$ vertices.\nLet $x$ be the $p$-dimensional random vector where each entry corresponds to a vertex.\nWe consider the following model that $x$ satisfies \n$$\nx_j = \\sum_{k\\neq j} \\beta_{jk}(z) x_k + \\epsilon_j\n$$\nwhere $\\beta_{jk}(z) \\neq 0$ if there is an edge between the vertex $j$ and $k$, $\\epsilon_j$ is an independent random variable whose expectation is $0$ and $z$ is a latent variable.\nWe further assume that $\\beta_{jk}(z)$ is parametrized by a neural network and the neural network is parametrized by $\\theta$.\nNamely, we rewrite $\\beta_{jk}(z)$ to be $\\beta_{jk}(z,\\theta)$.\n\nNow, suppose we are given $n$ samples $\\{(x^i,z^i)\\}$ generated by the aforementioned model.\nThe authors designed an algorithm to approximate the neural network parameter $\\theta$ and hence $\\beta_{jk}$.\nThe main result of this paper is to show that the generalization error decays in the rate of $\\frac{1}{n}$ (up to log factors).\nThe authors also provide the experimental results.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": ".", "broader_impact_concerns": "The work is mainly a theoretical result.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "FTpCLNQxKl", "reviewer": "Reviewer_FDHz", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- Extending beyond the Gaussianity assumption is an interesting and important topic in graphical model estimation.\n- The estimation algorithm is reasonable.\n- Finite samples bounds are provided for the method.\n- The empirical results are convincing.\n\nMain comments and/or weaknesses:\n- The function class considered seems to be rather restricted.\n- What is the benefit/motivation of assuming $f_{jk}(z,x_k)=\\beta_{jk}(z)x_k$? Is it precisely about improving statistical efficiency? Can't one just use the general method provided in https://www.jmlr.org/papers/v25/21-0022.html and https://arxiv.org/abs/2305.11379 to handle the more general case of $f_{jk}(z,x_k)$? This should be discussed to make clear the motivation of assuming a more restrictive functional class.\n- To make the motivation stronger, it would be better to explain more real-world examples about how the assumed function class will hold, in addition to the example of \"age\" in the introduction.\n- In p. 5, Is the assumption that $x$ and $z$ are bounded in $[-C,C]$ a common assumption? Doesn't that exclude the case where some of the variables are, e.g., Gaussian?\n- The rate/bounds in the theorems translate into Euclidean/Frobenius norm (e.g. as mentioned in Remark 6). How does such Frobenius norm consistency translate into bound on graph recovery which often demands more (e.g., see https://proceedings.neurips.cc/paper_files/paper/2008/file/61f2585b0ebcf1f532c4d1ec9a7d51aa-Paper.pdf for the Gaussian case)? If it cannot be translated, then the claim \"we establish guarantees on the recovered graph\" in the second contribution may not be completely accurate.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper aims to learn undirected graphical models where the structure is covariate-dependent, and develops an estimation method involving neural network. The method assumes specific function class where the linear coefficients depend (nonlinearly) on some covariate. The authors then provide theoretical guarantee for the finite-sample bounds of the method. Empirical results on synthetic and real data are provided.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See above section.", "broader_impact_concerns": "There are no broader impact concerns as far as I am concerned.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The authors have addressed all reviewers' concerns.
Decision: Accept (poster)