Bilevel Network Learning via Hierarchically Structured Sparsity
Authors: Jiayi Fan
Jingyuan Yang
Shuangge Ma
Mengyun Wu
Venue: NeurIPS 2025 poster ()
Keywords: neural network
network estimation
hierarchical selection layers
sparse regularization
Primary area: machine_learning_for_sciences

Abstract: Accurate network estimation serves as the cornerstone for understanding complex systems across scientific domains, from decoding gene regulatory networks in systems biology to identifying social relationship patterns in computational sociology. Modern applications demand methods that simultaneously address two critical challenges: capturing nonlinear dependencies between variables and reconstructing inherent hierarchical structures where higher-level entities coordinate lower-level components (e.g., functional pathways organizing gene clusters). Traditional Gaussian graphical models fundamentally fail in these aspects due to their restrictive linear assumptions and flat network representations. We propose NNBLNet, a neural network-based learning framework for bi-level network inference. The core innovation lies in hierarchical selection layers that enforce structural consistency between high-level  coordinator groups and their constituent low-level connections via adaptive sparsity constraints. This architecture is integrated with a compositional neural network architecture that learn cross-level association patterns through constrained nonlinear transformations, explicitly preserving hierarchical dependencies while overcoming the representational limitations of linear methods. Crucially, we establish formal theoretical guarantees for the consistent recovery of both high-level connections and their internal low-level structures under general statistical regimes. Extensive validation demonstrates NNBLNet's effectiveness across synthetic and real-world scenarios, achieving superior F1 scores compared to competitive methods and particularly beneficial for complex systems analysis through its interpretable bi-level structure discovery.

Bilevel Network Learning via Hierarchically Structured Sparsity
Jiyi Fan'; Jingyuan Shuangge Ma², Mengyun Wult
'Scoool of Statisties and Data Science, Shanghai University of Finance and Economics *Department of Biostatistics, Yale School of Public Health fanj sanjiayisstun sufe edu. en, jang.jiigguandeee shuangge ma@yall. edu, wu nu.mengyun@mai il shufe. thufe.edee edu. cn
sufe edu. en,
Abstract
Accurate network estimation serves as the cornerstone for understanding complex systems across scientific domains, from decoding gene regulatory networks in sys- tems biology to identifying social relationship patterns computational sociology. Modern applications demand methods that simultaneously address two critical challenges: capturing nonlinear dependencies between variables and reconstructing inherent hierarchical structures where higher-level entities coordinate lower-level components (e.g.. functional pathways organizing gene clusters). Traditional Gaus- sian graphical models fundamentally fail in these aspects due to their restrictive linear assumptions and flat network representations. We propose NNBLNet, a neural network-based learning framework for bi-level network inference. The core innovation lies in hierarchical selection layers that enforce structural con- sistency between high-level coordinator groups and their constituent low-level connections via adaptive sparsity constraints. This architecture integrated with a compositional neural network architecture that learn cross-level association patterns through constrained nonlinear transformations, explicitly preserving hierarchical dependencies while overcoming the representational limitations of linear methods. Crucially, we establish formal theoretical guarantees for the consistent recovery of both high-level connections and their internal low-level structures under general statistical regimes. Extensive validation demonstrates NNBLNet's N effectiveness across synthetic and real-world scenarios, achieving superior scores compared to competitive methods and particularly beneficial for complex systems analysis through its interpretable bi-level structure discovery.
Introduction
Network estimation is fundamental task across many disciplines such genetics, finance, and social science-whine uncovering the structure of dependencies among variables can yield critical insights into the underlying mechanisms of complex systems B 38 441 For instance, in genetics, gene regulatory networks can be reconstructed from analyzing multi-patient omics datasets to identify statistically significant edges between genes. In social science, co-authorship networks connect researchers through jointly published papers, with edge weights reflecting collaboration intensity, and the resulting network reveals meaningful patterns of scholarly connection and knowledge flow.
key characteristic shared across these domains is that variables can often be naturally organized into groups. Such group information is readily available as direct labels in practice: for instance, most genes have well-annotated pathway information in databases like KEGG, which naturally defines
"hese authors contributed equally.
Conresoonding author. Email: wu mengyundmaill shufe. edu. cn
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

functional groupings, while public social datasets commonly include demographic or professional categories such as gender or research field that provide meaningful group divisions. Under such a grouped organization, connections may occur both within these groups groups-rrfleecing functional or structural coherence and between
higher-order dependencies or cross-functional regulation (26 6]. Crucially, these systems exhibit hierarchical dependency architectures: depen- dencies among low-level variables often emerge only when their parent groups share systemic interdependencies. For instance, in genomics, genes within pathway rarely interact with genes in unrelated pathways unless cross-pathway regulatory mechanisms exist [4. TQ]. This hierarchical structure poses significant modeling challenges, conventional network estimation methods that ig- nore group-level level dependencies often fail to capture emergent system behaviors, yielding fragmentary results that lack systematic interpretability. This motivates bi-level network estimation, which aims to recover hierarchically structured dependencies at both the group and variable levels.
Gaussian Graphical Models (GGMs) are widely used for learning network structures by estimating the inverse covariance matrix under multivariate Gaussian assumptions due to their ability to capture conditional dependencies between variables
[39] Furthermore, a series of methods reformulate the GGM estimation problem into a set of sparse linear regression models, often offering greater computational efficiency (44). Extensions to bi-level GGMs have also been proposed to capture grouped variable associations [29]. However, GGMs are restricted to modeling linear relationships and rely heavily Gaussianity, limiting their applicability in real-world data with nonlinear or non- Gaussiusssa structure. To address this gap, a number of model-free or nonparametric approaches have been proposed [33 36 32] While these methods alleviate distributional assumptions, they often struggle to model highly complex interdependencies modern high-dimensional systems, where traditional kernel-based graphical techniques lack the representational capacity to capture intricate hierarchical patterns. Neural network (NN) has demonstrated superior performance in capturing nonlinear relationships compared to conventional nonparametric methods [18] Nevertheless, neural network-based bi-level network estimation remains unexplored.
We propose neural network-based framework for bi-level network estimation to address the chal- lenge of modeling complex, hierarchical dependencies in grouped variables. This work pioneers the use of neural networks for bi-level network inference, effectively capturing nonlinear interdepen- dencies both within and across groups, beyond the capacity of traditional linear models. The key contributions of our work are:
Bi-level Network Estimation: We propose a structured estimation framework that recovers hierarchical network architectures through dual-layer selection mechanisms, capturing bi- level dependencies between variables by incorporating group information and identifying group dependencies simultaneously.
Hierarchical Nonlinear Architecture: We integrate compositional neural network architec- tures to model cross-level dependencies via constrained nonlinear transformations, explicitly preserving hierarchical structures while overcoming the representational constraints of linear methods.
Theoretical Guarantees: We formally introduce the notion of bi-level selection consis- tency and establish a rigorous selection consistency result under high- nigh-dimensional regimes, ensuring the method to reliably recover the true network structure at both levels.
Empirical Validation: We validate our method through comprehensive experiments on both synthetic and real-world datasets. The results highlight the advantages of our approach in accurately estimating complex network structures. especially in scenarios where traditional GGM-based methods fail due to nonlinearity or distributional misspecification.
1.1 Notation
Consider zi independent observations with p variables partitioned into L predefined groups (Fig- ure A), where the 1-th group contains pr variables such that E Pt = p. For the i-th obser- The group mem- bership of the j-th variable is specified by C2 E (1....,4 L}. We further define two derived Ti Tiip)) representing the i-th observation with
vation, define the continuous measurement vector as x; (I.1....,
quantities: æi,-j = the j-th variable excluded and T. denoting the subvector containing variables from group | in the i-th observation. We define the operator @ as a block-wise product between vector and

---

variables
each C =L,
Zig = teu
Input
=hy(0: o
) +ay
Output
Zij ~ Zie Low-leve (variable)
High (group) alges
Group
Group
A. Data matrix
Node-wise regression
For each E (...p):: la (m.,)) = 77th 1o o.. o (r, (
L groups
Hierarchical network estimation:
High level: ou of
B.g
Low level: By r) 70 2tj 2ik
===
Pig
(a,()))))
hidden layers
Bi-level selection layers
D. NNBLNet structure:
(z)..))) estimated fa, (wi..)
Figure 1: Workflow of the proposed bi-level network estimation. (A) Input = (1I,1, ... Lipp) 's: observations with p variables partitioned into L groups. (B) Node-wise nonlinear regression for modeling the relationships between the jth variable Ti and the others -;,-j, as well as those between groups Lt, and ãz,r" Here, C₁ indicates that the jth variable belongs the Ith group. (C) Output: Bi-level network, including both low-level edges and high-level edges. (D) Network structure of the proposed NNBLNet, including the two hierarchical selection layers and multiple hidden layers for accommodating the nonlinear relationships among both low-level variables and high-level groups
partitioned vector (or matrix). Specifically, let A (aI..... am)
E
R mmxx be vector and
B = ((B")
B (m)
E RRaN be a matrix that can be decomposed into m con-
tiguous blocks, i.e., ((!) RRIXN for I l....,m with It= ni The block-wise product
A O B E R is defined as A o B
B
a2
B [
.....
B (m)
That is, each element aj scales all entries in the corresponding block (())
2 Neural Network-based Bi-level Network Estimation (NNBLNet)
Traditional methods, such as GGMs, typically handle linear conditional dependencies captured via precision matrices. However, realistic systems often exhibit complex, nonlinear associations that linear models may not adequately capture. To overcome these limitations, our framework formalizes nonlinear relationships among variables through node-wise nonlinear regression. Specifically, for j-th variable with group assignment C = consider
Fij = (x Eij h, ( (STT....,
T.
Eij.
where T..-, specifically excludes the j-th variable from its native group I, Eij is independent and identically distributed sub-Gaussian stochastic noise with mean zero and sub-Gaussian parameter

---

and hy RP-1 + is sparse nonparametric mapping (Figure 1 B). Model (1) captures the nonlinear relationships between zi,j and all other variables z; enabling more accurate characterization of complex system mechanisms than linear approximations permit. The group information is incorporated in (I) through structured group inputs øtg. This strategy can also accommodate nonlinear high-level relationships between group / and all other groups, leading to deeper investigation of the underlying structures (Figure
For estimating h-)() and performing bi-level network estimation, we propose a neural network with two hierarchical selection layers and D hidden layers (Figure D):
Ja (2), T. woTT o.. og. (r, (0, 8 3 F
(2)
with aj { [1 (d) .b |d=1] being the parameter Here, each affine transformation given by (1)(((()) (u) و (d) + 050) with learnable parameters A (d)
E IR x and b,²)
R=)"xx where w is the width of the d-th hidden layer. The activation function @ is chosen as the ReLU function and the block-wise product defined in Notation.
In 2 we innovatively introduce two selection layers. Specifically, the first selection layer involves the vector 0, = (011.... ........... Ou...., 0nl) which encodes group-level (high-level) dependencies with
0 = 1. Simultaneously, the second selection layer involves the matrix T3 R Re x(p-1), which captures specific associations between variable j and the other p variables. Hence, by integrating these hierarchical components, the relationships among low-level variables are explicitly characterized by 0; @ r where high-level group dependencies act as latent scaffolds that constrain enable low-level dependencies.
(1) ( 1" Denote a = ( {
12=1 (,,
with O = 02,, 0r), and r the k-th
1
column of j. Integrating Sa, (1i,-1), .... Sa, (2))-p), we propose a dual-penalized estimator designed to achieve bi-level sparsity for network estimation:
On = argmin
Σ
I(aj, Fz.j. Di, xi,- + A
Σ t<P
10
+ A2
(3)
J=1
with ((ay, Fiiy Li, = (xig - (xx, and ||. being the I--norm of a vector.
Here, the first term constitutes quadratic reconstruction loss. The second term implements a Lasso penalty on 0g, inducing sparsity through element-wise shrinkage towards zero, which enables automatic high-level edge selection. This foundational assumption of inter-group sparsity aligns with established practices in prior methodological work [5] 29 This is well illustrated the context of
genomic networks, where pathway-to huthway-to-patwayy connections exhibit natural sparsity since regulatory relationships occur only between specific pathway pairs rather than universally [27]. The third
term employs group Lasso penalty on !! enforcing simultancous shrinkage of entire parameter vectors to zero, thereby facilitating low-level edge selection. The tuning parameters a1 and X govern edge selection stringency, with larger values yielding sparser network structures.
Based on o., to mitigate regularization bias and enhance selection consistency, we further introduce an adaptive bi-level sparse estimator:
M
10 10
on = arggin
(a Fig Lig. + +1
+
C2
[.
T:
Ier
j=1
with Çi and (2 being two tuning parameters for controlling network sparsity and y being a positive constant. Here, the weights |/ and [.A] connections with larger initial estimates while amplifying the shrinkage for smaller ones, resulting in more accurate edge selection owing to reduced estimation bias.
are introduced to adaptively reduce penalties for
The final bi-level network is constructed based on ã.. Specifically, high-level connections between groups | and L' are established when either pur + 0 or dvi * 0. This high-level connectivity then

---

informs low-level interactions: an edge forms between node j in group I and node k in group I' if either Ou / + or f extended to accommodate overlapping groups, with the details provided in the Appendix.
& + O. The proposed bi-level network estimation strategy can also be
3 Statistical Properties
3.1 Approximation Error
Without loss of generality, assume that w, E X [0, [J. Define the 3-Hölder smooth class:
H. (10, 11" Bo)
f: [0. R. max 1 flloo S Bo, max sup j0" f(x)-0" 0000-00fp00 f(y)l [|]
I Bo
|| zÿμ 11e-vll (2)
where 0 = 80r .... -0. with a
"u.,...,a,)] E NJ and No denotes the set of non-negative integers, ||a||1 = Ei=1 |ag| and [8] denotes the largest integer strictly smaller than 3. If function Bo),
belongs to ([0, Bo), then all the partial derivatives up order [6] exist, and the partial derivatives of order [ are B - [8] Hölder continuous.
Asummption 3.1. For j = ,,. the target sparse function h, resides in the sy-sparse function class Faj, where F, = {h: [0, 1jm-1 - [0, 3h: [0, [0, 1] E ([0, 1]", Bo) s.l. h(x)) = h(t)), Va E [0, 1]1-1]] with ă E IR representing the relevant s-dimensional subvector of x.
Assumption 3.1 posits that the target functions exhibit sparsity and smoothness. This is common requirement neural network approximation theory [15 4 In addition, the assumption is grounded in reality: for instance, in gene regulatory networks, genes are sparsely rather than fully connected, and nonlinear effects vary smoothly rather than abruptly.
Next, we first establish the approximation capabilities of our proposed sparse ReLU feedforward neural network architecture for sparse nonlinear functions. We define the width of neural network as the maximum width among its hidden layers.
Theorem 3.2 (Approximation Error). If h(x)) FFrr then for any positive integer N and M, there exists sparse ReLU feedforward neural network f(x) with width W = 38(13] 1)23ss |ا8 N (logg(8N)) depth D = 21([3] + 1)2 log,(8M) + 2s, and the selection layer parameter & ) being able to rearranged as [u (u=") xs x(p-1=x)))) such that:
1. Support equivalence: both h(x)) and f(x) are s-sparse functions and have the same support:
2. Uniform approximation:
sup
[f(x) h(x)) < (1000(0))1+))000000000)0((0)1110
(5)
|| |
Under Assumption Theorem 3.2 indicates that sparser nonlinear functions require smaller network dimensions for accurate approximation an observation aligning with intuition. Therefore, the target sparse function hy can be approximated by sparse ReI.U feedforward network with parameter set aj, denoted as faz.
3.2 Bi-level Selection Consistency
Let cy NN - be the column label mapping where (k) denotes the column T2 corresponding to the k-th variable, that is, ()(a) = & if k < and (y(k) = h 1 if k > j. Denote [p] and [Z.] as the set {1,... .P} and (1, . L}. respectively, and A as the complement of A. In
addition, we formalize three fundamental sets: A {A [p] 0 0 (gl:s,(())) + +0).
By k E (p) r (+))) = 0). and P: l' [L] 8 = =).
Then, we can define the bi-level selection consistency as follows.
Definition 3.3 (Bi-level Selection Consistency). We say the estimator (0, (T,})=1) achieves bi-level selection consistency v8 > 0, 3N3 IN such that Vn > Ng:
5

---

mental equivalence: 1. High-level true positives: Bw * foror all (1,1) [L] Pf:
2. High-level true negatives: Bu foror all (1,1) E [4] Pi;
(oc,orr,)
3. Low-level true positives:
* for all (j,k) [p] x Aj;
(0cc, Be((K)
4. Low-level true negatives:
for all (j.h) E (p) x A3
with probability at least 1 5.
Let a consist of all aj, o and W denote the feasible parameter space with a' e W. We define the population risk R(a)) and empirical risk (ax associated with the squared error as:
R(a) = Σ Rg(aa) j=1
Σ j=l
E (2) -x,))² =²1
(6)
Ra(a) Σ Rnj(as) 1
(Ja, (2i -za))
1=1
The optimal parameter is defined as:
H fa E W: R(a)) R(a"))
Despite the intricate geometric structure of H(" we establish in Lemma A.1 (Appendix) the funda-
0 E E
Saa, = Sag: Vj (1,.....
This equivalence implies that both the estimation and selection consistency can be simultancously attained by controlling the proximity of the parameters to the optimal set H.. Define d(a,, H²) inføeи- ||a A||-
To establish these theoretical guarantees, we introduce the following assumptions.
Asummption 3.4. Define as a class of ReLU feedforward neural networks f& 10,1 [0, with parameter o depth D, width W, size S (the number of elements in &e) and B- B-Lipschitz continuity. We assume that fa, E F for all j 1.2......
Asummption 3.5. There exist C9 > and v >> 2 such that R(33) - -R(x*)) 2 crd(B. y for all ( W.
Assumption 3.4 posits the boundedness of the ReLU neural network function class. This assumption holds in practice, real-world data (e.g.. gene expression levels) and weights are naturally bounded by physical constraints. Assumption 3.5 is a technical assumption. For fixed network with an analytic activation function, it holds and can be justified by Lojasiewicz's inequality 99.
Under Assumptions 3.4 and 3.5 the following convergence properties hold:
Theorem 3.6 (Group Lasso Lasso Convergence). Let p o(log n), SDlog(S) O(n O(n-h) b), and X2 On 3), then there exist c3 > 0 and C₁ such that
log n n
d(a,, H²) < cg
(10)
and
|
(R))
Σ
10 s logn
log n n
(11)
KEB,
1=1 Fire
holds with probability at least 81 with 01 4n (nt +
log
(32ent)
(logn)nt
log

---

Theorem 3.7 (Consisteney of the adaptive bi-level sparse estimator). Let > 0, > 0, = O(n )) with adaptive bi-level sparse penalty has bi-level selection consistency, and there exists C2 > Osuch that
22 = O( ) 0 15 15*)) and E2 = O (n कu 15=t then the estimator Ox.
d( S csn (-
(12)
holds with probability al least - 81-
Under Assumption 3.4 and Theorem 3.7 the estimation error can be bounded by leveraging the Lipschitz continuity Ifa, faz S Begn "y This result implies that the estimation error decays as sample size n increases. Combined Theorem 3.2 and Theorem 3.7] we can easily establish
the estimation consisteney of fa, by triangle inequality. That is, [fa, hyl S Boyn wo- (- te)/v 19 [[3] ())3++11+(())/)))((NNN))) -M8//// 228/* with probability least 82.
Computation
The proximal gradient descent algorithm for training NNBLNet is outlined in Algorithm (see Appendix). Following the convergence conditions in Theorems 3.6 and 3.7 we set X = S
e.n e-n-1/8 for = 1,2, with e as a tunable constant. This choice follows conventional rate settings in high-dimensional sparse inference (22]. Empirically, satisfactory accuracy across diverse data settings. For finer calibration, these parameters can be selected via cross-validation (see Appendix|A.A. A.4). Consistent with common practice, we set y 1. Regarding neural network hyperparameters. exploratory experiments (Appendix [A.4) showed that configuration of 1000 training epochs, three hidden layers, and 50 nodes per layer offers an optimal balance between accuracy and computational efficiency. This fixed setup was used in all experiments to ensure reproducibility and consistency, in line with common practice 33
recommend
0.35, as it achieves
Experiment
To comprehensively validate our methodology, established a dual evaluation framework compassing both synthetic benchmarks and real-world networks. The synthetic analysis included systematically constructed networks with predefined nonlinear and linear patterns, enabling con- trolled assessment of relationship modeling capabilities. For real-world validation, we analyzed four distinct network types: (I) Friendship- social connections among high school students from a high school; (2) Co-authoshipp-oollaoraii relationships in academic publications; (3-4) BRCA and LUAD -ene regulatory networks derived from The Cancer Genome Atlas (TCGA) breast cancer and lung adenocarcinoma data, respectively.
In addition the proposed NNBLNet, we conducted systematic comparisons with three alternative methodologies: (1) BGSL 7. Bayesian Gaussian graphical modeling framework with explicit group-structured variable representations; (2) Fair Glasso (25). which enhances sparse precision matrix estimation through group-aware regularization constraints; and (3) DeepGRNCS [19]. a multi-task deep learning architecture specializing in joint network inference. Experimental details, including the synthetic settings and introductions to the real-world datasets, are provided in Appendix A.5. To supplement the results presented in the main text, we also conducted a series of sensitivity analyses in Appendix A.6] including an assessment of model generalization capability based on the stability of network estimation and performance comparisons under different sample sizes, group sizes, and group label misclassification rates, as well as downstream analysis of the LUAD dataset.
5.1 Synthetic Networks with Nonlinear and Linear signals
In [ NNBLNet consistently outperformed all baseline methods in nonlinear scenario, achieving the highest F1-score of 0.772. This improvement was attributed to its hierarchical design, which jointly modeled both within-group and between-group edges through bi-level regularization. Fair Glasso and BGSL, which impose sparsity within the GGM framework, yielded comparable results with F1 -cores of 0.710 and 0. 0.709, respectively. DeepGRNCS, a deep learning- method, performed better than classical approaches but fell short of our proposed method, suggesting that
7

---

integrating bi-level selection is crucial for recovering multiscale network structure. These results demonstrate the efficacy NNBLNet in uncovering meaningful networks with greater accuracy and structural coherence.
Under the linear setting, Fair Glasso attained the highest Fi-score (0.790), aligning with its design premise for sparse Gaussian graphical models. NNBLNet remained highly competitive (FI = 0.769), which illustrates its adaptability across both linear and nonlinear data regimes. Although specifically designed for hierarchical nonlinear dependencies, it maintained robust performance even in simpler settings. DeepGRNCS also delivered favorable results, while BGSL achieved a balanced recall- precision trade-off. These findings verify that NNBLNet not only performs stably in conventional linear estimation tasks but also offers flexibility more complex environments.
Table Performance comparison of different methods for the synthetic network with nonlinear and linear relationships: Mean (SD) over 100) replicates.
Patterns
Metric
NNBLNet
Fair Glasso
BGSL
DeepGRNCS
Recall
0.872(0.016)
0.779(0.014)
0.790(0.017)
0.809(0.016)
Nonlinear
Precision
0.693(0.022)
0.653(0.021)
0.644(0.018)
0.681(0.018)
Fi-score
0.772(0.014)
0.710(0.017)
0.709(0.013)
0.731(0.017)
Recall
0.881(0.018)
0.901(0.016)
0.836(0.016)
0.846(0.016)
Linear
Precision
0.675(0.020)
0.704(0.019)
0.692(0.019)
0.671(0.020)
Fl-score 0.765(0.013)
0.790(0.013)
0.757(0.015)
0.763(0.017)
5.2 Four Real-World Networks
We further evaluated our method on four real-world datasets spanning social networks and biological systems: Friendship, Co-authorship, BRCA, and LUAD.. Performance results summarized in T
Table 2: Performance comparison of different methods for the four real-world networks.
Dataset
Metric
NNBLNet
Fair Glasso
BGSL.
DeepGRNCS
Recall
0.875
0.859
0.798
0.811
Friendship
Precision
0.804
0.703
0.735
0.686
Fl
0.838
0.771
0.765
0.745
Recall
0.712
0.643
0.622
0.655
Co-authorship
Precision
0.678
0.606
0.591
0.610
Fl
0.695
0.619
0.604
0.632
Recall
0.764
0.812
0.676
0.742
BRCA
Precision
0.618
0.652
0.552
0.601
F1
0.683
0.723
0.607
0.664
Recall
0.641
0.625
0.597
0.654
LUAD
Precision
0.526
0.517
0.489
0.539
Fl
0.577
0.566
0.537
0.591
Across all four datasets, NNBLNet consistently achieved either the best or second-best Fl-score. demonstrating its versatility and robustness in diverse application domains. On the Friendship network, it significantly outperformed competitors with an Fl-score 0.838, reflecting its strength in modeling community-driver structures. In the Co-authorship setting, where overlapping communities and latent hierarchies were expected, NNBLNet again led with an FI-score 0.695, significantly outperforming Fair Glasso (0.619), BGSL (0.604), and DeepGRNCS (0.632). These results indicate the advantage jointly modeling inter-group and intra-group dependencies.
For the biological datasets, BRCA and LUAD, NNBLNet maintained top-tier performance. While Fair Glasso achieved a slightly higher F1 score in BRCA (0.723), NNBLNet offered more balanced precision-recall trade-offs, especially in LUAD, where it achieved the second highest FI-scon

---

(0.577). The improvement in biological contexts suggests that our method effectively captures complex relationships between molecular factors, including both nonlinear individual factor-level dynamics and pathway-level crosstalk.
Collectively, these results demonstrate that NNBLNet is highly competitive and generalizable frame- work for network estimation, performing consistently well across domains with varying complexity and noise characteristics.
5.3 Ablation Study
To validate the necessity of the proposed bi-level architecture, we performed ablation studies com- paring our model against two simplified variants: the flat-structured NNNet, which ignores the bi-level design, and Modified-NNNet. which uses group labels but omits intra-group adjacency modeling. As shown in Table NNBLNet consistently outperformed both variants on synthetic and real-world datasets. The significant performance gain confirms that the bi-level structure essential for achieving higher recall, precision, and FI scores.
To further isolate the benefits of the bi-level architecture beyond group information alone, we conducted sensitivity analyses in scenarios with extreme sparsity of inter-group edges (Table 4 where the ratio of inter-group edges to total edges (7) was set to 0%, 0.5%, and I %. The results showed that NNBLN Nett consistently achieved the highest F1 by balancing recall and precision. In addition, the modified-NNNF -NNNet failed to capture inter-group edges and exhibited declining recall as increased, while the NNNet, which ignored grou- sparsity, produced excessive false positives under sparse settings. These findings confirmed that the bi-level hierarchical sparsity structure is essential for accurate network inference.
Table 3: Performance comparison of NNBLNet, NNNet, and Modified-NNNet across six datasets (mean values and standard deviation for synthetic datasets).
Dataset
Method
Recall
Precision
F1
Nonlinear
NNBLNet
0.872 (0.016)
0.693 (0.022)
0.772 (0.014)
NNNet
0.844 (0.017)
0.656 (0.019)
0.738 (0.014)
Modified-NNNNt
0.782 (0.020)
0.660 (0.021)
0.716 (0.017)
Linear
NNBLNet
0.881 (0.018)
0.675 (0.020)
0.769 (0.013)
NNNct
0.857 (0.014)
0.628 (0.017)
0.727 (0.013)
Modified-NNNet
0.802 (0.018)
0.627 (0.019)
0.704 (0.016)
Friendship
NNBLNet
0.875
0.804
0.838
NNNet
0.842
0.745
0,790
Modified-NNNe
0.774
0.734
0.753
Co-authorship
NNBLNet
0.712
0.678
0.695
NNNet
0.674
0.637
0.649
Modified-NNNa
0.612
0.636
0.624
BRCA
NNBLNet
0.764
0.618
0.683
NNNet
0.709
0.563
0.627
Modified-NNNG
0.655
0.566
0.607
LUAD
NNBLNet
0.641
0.526
0.577
NNNet
0,613
0.502
0.552
Modified-NNNB
0.552
0.503
0.527
Discussion
This work presents NNBLNet, neural network framework for bi-level network inference. The method is built upon a key structural prior that represents dependency structures complex systems as hierarchical. This hierarchical modeling is implemented through two mechanisms: intra-group information sharing, which amplifies weak signals via latent pooling within groups, and inter-group sparse transmission, which gates cross-group connections through switches (0w) to suppress imree- vant noise. NNBLNet represents paradigm shift structured network inference by unifying neural network representation learning with hicrarchical sparsity constraints. Its innovative hierarchical

---

Table 4: Performance of NNBLNet, NNNet, and Modified-NNNe. under sparse grour edges (mean and standard deviation over 100 replicates).
Setting
1
Method
Recall
Precision
FI
Nonlinear
0%
NNBLNet
0.752 (0.018)
0.684 (0.020)
0.716 (0.017)
NNNet
0.744 (0.019)
0.619 (0.023)
0.676 (0.018)
Modificd-NNNet
0.737 (0.018)
0.682 (0.020)
0.708 (0.016)
0.5%
NeBLNet
0.760 (0.017)
0.696 (0.019)
0.727 (0.016)
NNNet
0.752 (0.018)
0.616 (0.022)
0.677 (0.017)
Modificd-NNNct
0.717 (0.019)
0.673 (0.021)
0.694 (0.017)
1%
NNBLNet
0.774 (0.016)
0.715 (0.018)
0.743 (0.014)
NNNet
0,764 (0.018)
0.615 (0.023)
0.681 (0.019)
Modificd-NNNL NNNNet
0.696 (0.020)
0.663 (0.022)
0.679 (0.019)
Linear
0%
NNBLNet
0.740 (0.019)
0.672 (0.021)
0.704 (0.018)
NNNet
0.732 (0.020)
0.607 (0.024)
0.664 (0.019)
Modified-NNNG.
0.724 (0.019)
0.671 (0.021)
0.696 (0.017)
0.5%
NNBLNet
0.750 (0.018)
0.684 (0.020)
0.715 (0.017)
NNNet
0.741 (0.019)
0.605 (0.023)
0.666 (0.018)
Modified-NNNee
0,706 (0.020)
0.662 (0.021)
0.682 (0.018)
1%
NeBLNet
0.766 (0.017)
0.702 (0.019)
0.733 (0.015)
NNNet
0,753 (0.019)
0.605 (0.024)
0.671 (0.019)
Modified-NNNet
0.681 (0.021)
0.652 (0.022)
0.666 (0.020)
selection layer explicitly captures bi-level dependencies: local associations between individual vari- ables and global coordination among groups. This architecture effectively models complex nonlinear relationships while maintaining interpretability. Theoretically, we establish a bridge between neural modeling and statistical guarantees by proving estimation consistency and exact bi-level selection consistency. Empirical evaluations and real-data analyses demonstrate NNBLNet's ' effectiveness, showing it achieves superior Fl scores compared competing methods.
6.1 Limitation
This study has several limitations that point to valuable directions for future research. Theoretically, our analysis relies on the assumption of sub-Gaussian noise, which, while common in statistical literature, is often difficult to verify in practical applications. Despite this limitation, empirical comparisons against authoritative ground-truth networks demonstrate that our method still achieves competitive network reconstruction accuracy, thereby offering partial validation of its practical effectiveness. Moreover, our method relies on pre-defined group labels, and performance may decline if these noisy incomplete. Future extensions could jointly infer group memberships and associations or incorporate unsupervised techniques (c.g., spectral clustering) when prior labels are unavailable, enhancing robustness and applicability across diverse scenarios.
6.2 Broader Impact
NNBLNet bridges neural network learning and hierarchical sparse inference to enable interpretable bi-level network discovery across scientific domains, from identifying gene-pathway interplay in disease mechanisms modeling individual-grous dynamics in social systems. Its adaptive sparsity design balances predictive power with mechanistic interpretability, offering actionable insights for precision medicine and policy-making.
Acknowledgments and Disclosure of Funding
We thank the Area Chair and the anonymous reviewers for their insightful feedback, which was instrumental in strengthening this paper. This research was supported by the National Natural Science Foundation of China (12071273); MOE Project of Humanities and Social Sciences (25YJCZH291); Shanghai Rising-Star Program (2201403500); Shanghai Science and Technology Development
10

---

Funds (23JC1402100); Shanghai Research Center for Data Science and Decision Technology; National Institutes of Health (CA204120); and National Science Foundation (2209685).
References
[1] M. Anthony and I.. Bartlett. Neural network learning: Theoretical fleedratial Cambridge University Press, 2009.
[2] P. L.. Bartlett, N. Harvey, C. Liaw, and A. Mchrabian. Nearly-tight ve-dimension and pseudodi- mension bounds for piecewise linear neural networks. Journal of Machine Learning Research, 20(63):1-17, 2019.
[3] Chen, L. M. Ashwood, O. Kondrashova, A. Strasser, G. Kelly, and K. D. Sutherland, Breathing new insights into the role of mutant p53 in lung cancer. Oncogene, 44(3):115-129) 2025.
[4] Y. Chen, Q. Gao, F. Liang, and X. Wang. Nonlinear variable selection via deep neural networks. Journal of Computational and Graphical Statistics, 30(2):484 492, 2021.
[5] L. Cheng, L. Shan, and I. Kim. Multilevel gaussian graphical model for multilevel networks. Journal of Statistical Planning and Inference, 190:1-14,2011
[6] I. Choi and W. C. Kim. A multifaceted graph-wise network analysis of sector- based financial instruments' price-based discrepancies with diverse statistical interdependencies. The North American Journal of Economics and Finance, 75:102316, 2025.
[7] A. Colombi, R. Argiento, L. Paci, and A. Pini. Learning block structured graphs in gaussian graphical models. Journal of Computational and Graphical Statistics, 33(1):152-165, 2024.
[8] L. Cowen, T. Ideker, B. J. Raphael, and R. Sharan. Network propagation: a universal amplifier of genetic associations. Nature Reviews Genetics, 18(9):551-562. 2017.
[9] V. C. Dinh and L. S. Ho. Consistent feature selection for analytic deep neural networks. Advances Neural Information Processing Systems, 33:2420-2433 2020.
[10] R. Feng. J. Gong. L. Wu, L. Wang, B. Zhang, G. Liang. H. Zheng, and H. Xiao. Mapk and hippo signaling pathways crosstalk via the raf- 1/mst22 interaction in malignant melanoma. Oncology Reports, 38(2): 189()):119912209) 1205, 2017.
[11] A. Halama, S. Zaghlool, G. Thareja, S. Kader, W. Al Muftah, M. Mook-Kanamori, H. Sarwath, Y. A. Mohamoud, N. Stephan, S. Ameling, al. A roadmap the molecular human linking multiomics with population traits and diabetes subtypes. Nature Communications, 15(1):1111, 2024.
[12] A. Herreros- Pootros, J. D. de Maya-Girones, S. Calabuig-Farifas, R. Lucas, A. Martínez, J. M. Pardo-Sánched Sánchez, S. Alonso, A. Blasco, R. Guijarro, M. Martorell, et al. Lung tumorspheres reveal cancer stem cell-like properties and a score with prognostic impact in resected non-small non-small-cell lung cancer. Cell Death & Disease, 10(9):660, 2019.
[13] V. A. Huynh-Thu, A. Irrthum, L. Wehenkel, and P. Geurts. Inferring regulatory networks from expression data using tree-based methods. PloS One, 5(9)e12777,, 2010.
[14] J.-W. Jang, M.-K. Kim, and S.-C. Bae. Reciprocal regulation of yap/taz by the hippo pathway and the small gtpase pathway. Small GTPases, 11())280--888 2020.
[15] Y. Jiao, G. Shen, Y. Lin, and J. Huang. Deep nonparametric regression on approximate manifolds: Nonasymptotic error bounds with polynomial prefactors. The Annals of Statistics, 51(2):691-716, 2023.
[16] M. Kanehisa, M. Furumichi, Y. Sato, Y. Matsuura, and M. Ishiguro-Watanabe. systems database as a model of the real world. Nucleic Acids Research, 53(D1))D672-D67, 2025.
Watanabe. Kegg: biological
11

---

[17] N. M. LaMarche, S. Hegde, M. D. Park, B. B. Maier, L. Troncoso, J. Le Berichel, P. Hamon, M. Belabed, R. Mattiuz, C. Hennequin, al. An il-4 signalling axis in bone marrow drives pro-tumorigenic myelopoiesis. Nature, 625(7993): 166-176,4417111 2024.
2024.
[18] Y. LeCun, Y. Bengio, and G. Hinton. Deep learning. Nature, 521(7553) 2015.
[19] Y. Lei, X.-T. Huang, X. Guo, K. Hang Katie Chan, and L. Gao. Deepgmec: deep learning- framework for jointly inferring gene regulatory networks across cell subpopulations. Briefings in Bioinformatics, bba4)33b 2024.
[20] A. Lewkowycz and G. Gur-Ari. On the training dynamics of deep networks with 1_2 regulariza- tion. Advances in Neural Information Processing Systems, 33:4790-4799, 2020.
[21] B. Liu, Y. Wei, Y. Zhang, and Q. Yang. Deep neural networks for high dimension, low sample size data, In IJCAI, volume 2017, pages 2287-2293, 2017.
[22] H. Liu and L.. Wang. TIGER: A tuning-insensitive approach for optimally estimating Gaussian graphical models. Electronic Journal of Statistics, 11(1):241 294, 2017.
[23] G. Mao, Pang, K. Zuo, Q. Wang, X. Pei, X. Chen, and J. Liu. Predicting gene regulatory links from single-cell ma-seq data using graph neural networks. Briefings in Bioinformatics, 24(6):bbad414, 14, 2023.
[24] N. Meinshausen and P. Bühlmann. High-dimensional graphs and variable selection with the Lasso. The Annals of Statistics, 34(3):1436 1462, 2006.
[25] M. Navarro, S. Rey, A. Buciulea, A. G. Marques, and S. Segarra. Fair glasso: Estimating fair graphical models with unbiased statistical behavior. Advances in Neural Information Processing Systems, 2024.
[26] M. A. Reyna, D. Haan, M. Paczkowska, L. P. Verbeke, M. Vazquez, A. Kahraman, S. Pulido- Tamayo, J. Barenboim, L. Wadi, P. Dhingra, et al. Pathway and network analysis of more than 2500 whole cancer genomes. Nature Communications, 11(1):729, 2020.
[27] Sanchez- M. Mina, J. Armenia, W. K. Chatila, A. Luna, K. C. La, S. Dimitriadoy, D. Liu, H. S. Kantheti, S. Saghafinia, et al. Oncogenic signaling pathways in the cancer genome atlas. Cell, 173(2):3213333 2018.
[28] L. Seninge, I. Anastopoulos. H. Ding. and J. Stuart. Vega is an interpretable generative model for inferring biological network activity in single-cell transcriptomics. Nature Communications, 12(1):5684, 2021.
[29] I.. Shan, 7. Qiao, L. Cheng, and I. Kim. Joint estimation of the two-level gaussian graphical models across multiple classes. Journal of Computational and Graphical Statistics, 29(3):562- 579, 2020.
[30] B. Shekar and G. Dagnew, LI Ll-regulated feature selection and classification of microarray cancer data using deep learning, In Proceedings of 3rd International Conference Computer Vision and Image Processing: CVIP 2018, Volume 2. pages 227-242. Springer, 2020.
[31] D. Stewart, M. Javadi, M. Chambers, C. Gunsolly, G. Gorski, and R. C. Borghaei. Interleukin- 4 inhibition of werrleukin---inducee expression of matrix metalloproteinase-3 (mmp-3) is independent of lipoxygenase and ppart activation in human gingival fibroblasts. BMC Molecular Biology, 8(1): 12, 2007.
[32] C. Su, Z. Xu, X. Shan, B. Cai, H. Zhao, and J. Zhang. Cell-tyee-speiffic co-expression inference from single cell ma-sequencing data, Nature Communications, 14(1):4846, 2023.
[33] T. Sun, Y. Wei, W. Chen, and Y. Ding. Gemome-wide association audy-based deep learning for survival prediction. Statistics in Medicine, 39(30) 2020.
[34] Y. Sun, Q Song, and F. Liang. Consistent sparse deep learning: Theory and computation. Journal of the American Statistical Association, 117(540):198111195 2022.
12

---

[35] D. Szklarczyk, K. Nastou, M. Koutrouli, R. Kirsch, F. Mehryary, R. Hachilif, D. Hu, M. E. Peluso, Q. Huang. T. Fang, et The string database in 2025: protein networks with direction- ality of regulation. Nucleic Acids Research, 5X(D1):D730-D7737, 2025.
[36] X. Wang, D. Choi, and K. Roeder. Constructing local cell-specific networks from single-cell data. Proceedings of the National Academy of Sciences, 118651))22113178118,, 2021.
[37] Y. Wang, T. Ran, Y. Li, L. Tian, L. Yang. Z. Liu, and B. Yao. Identification of jun gene and cellular microenvironment in response pd-1 blockade treatment in lung cancer patients via single-cell ma sequencing. Aging (Albany NY), 16(12): 1032), 2024.
[38] R. Xu, C. Mi, R. Mierzwiak, and R. Meng. Complex network construction of internet finance risk. Physica A: Statistical Mechanics and its Applications, 540:122930, 2020.
[39] M. Yuan and Y. Lin. Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19-35, 2007.
[40] M. Zhang, Y. Sun, and F. Liang. Sparse deep learning for time series data: theory and applications. Advances in Neural Information Processing Systems, 36:54509-54539, 2023.
[41] T. Zhang and H. Zou. Sparse precision matrix estimation via lasso penalized D-trace loss. Biometrika, 101(1): 1000 2014.
[42] Z. Zhang, H. Li, S. Jiang, R. Li, W. Li, H. Chen, and X. Bo survey and evaluation of web-based tools/databases for variant analysis of tega data. Briefings in Bioinformatics, 20(4):1524 1541, 2019.
[43] Y. Zhou, X. Wang, Li, M. Fan, M. Wang, W. Gao, Z. Shi, Wang, J. Tan, and F. Liu. positive feedback loop between serpinh1 and mmp-9/tf-f8 promotes lung adenocarcinoma progression. Cell Death Differentiation, pages 1-16, 2025.
[44] W. Zou, X. Hu, Z. Pan, Li, Y. Cai, and M. Liu. Exploring the relationship between social presence and learners' prestige in mooe discussion forums using automated content analysis and social network analysis. Computers in Human Behavior, 115:10 2021.
13

---

A Appendix Supplemennal Material
A.1 Related Works
In recent years, the explosive growth of high-dimensional data, such as high-throughput omics data, has generated unprecedented volumes of complex data, thereby elevating the importance of network estimation research. Among the various existing approaches, Gaussian Graphical Models (GGMs) have emerged as one of the most widely used techniques. In the GGM framework, network structures are inferred by estimating a sparse precision matrix (i.e., the inverse of the covariance matrix), which encodes conditional dependencies among variables conditional on all others [T] This estimation procedure can be formulated as series of sparse node-wise wise linear regressions. Compared to marginal or unconditional correlation-basee methods (e.g.. Pearson correlation), this conditional strategy offers more holistic view of system-level dependencies, potentially leading to more meaningful interpretations.
In real-world scenarios, networks often exhibit a bi-level hierarchical structure. This structure implies that some variables are organized into higher-level groups, with lower-level variables nested within these groups. For example, a gene pathway consists of multiple genes that collaborate to perform specific cellular or physical function. In this context, pathways represent the higher-level groups, while the individual genes within those pathways are the lower-level variables. It is also important to note that these groups are not independent of each other. To tackle hierarchical network estimation, several GGM- -bsed methods have been developed. Cheng et al. (2017) [5] introduced a multilevel Gaussian graphical model for nested data structures. Shan et al. (2020) 29 proposed a framework for joint estimation of two-level GGMs across multiple classes. Colombi et al. (2024) 7 focused on learning block-structured graphical models using variable groupings. Notably, Fair Glasso 25 specifically leverages group information to estimate graphical models with provably unbiased statistical behavior, addressing fairness concerns in network inference. However, these GGM based methods can only capture linear dependencies and may ignore complex nonlinear relationship in the real world.
Deep neural networks (DNNs) and related machine learning models have gained widespread attention due their strong capacity for nonlinear approximation and representation learning, particularly in high-dimensional settings (21. These models excel at uncovering complex associations within large datasets, making them particularly appealing for network inference tasks. With the increasing availability of large-scale high-dimensional data, deep learning has become a cornerstone for network estimation. Researchers have recently developed several deep learning-basse methods that aim to reconstruct latent networks. These methods harness the expressive power of neural networks to model intricate dependencies between variables, often using architectures such as pre-trained deep neural networks [19], Variational Autoencoders (VAEs) [28) and Graph Convolutional Networks (GCNs) [23]. While promising, these approaches face several important limitations. First, most supervised and require labeled data such as group-specific regulatory annotations curated databases which are expensive and time-consuming to acquire. Consequently, the size and diversity of training datasets remain limited. Moreover, deep learning methods often lack theoretical guarantees and interpretability, which hampers their adoption in sensitive scientific domains.
To improve interpretability, the research community has proposed a variety of strategies aimed at making deep learning models more transparent and reliable. Among these, regularization-ben approaches have shown particular promise. By incorporating additional constraints (e.g., LI or L2 penalties), regularization can limit model complexity, promote sparsity, and enhance feature selection 30 22). This is especially valuable for network inference, where identifying the most influential variables (e.g., hub nodes) is often a primary goal. In scientific applications, regularization enhances both statistical reliability and practical interpretability. Recent theoretical work further supports the effectiveness of regularization for consistent variable selection in deep models [34 40), providing rigorous foundation for interpretable deep learningbbase network estimation. These theoretical insights not only deepen our understanding of model behavior but also inform the development of robust and generalizable algorithms.
A.2 Extension to Overlapping Groups
In many scientific applications, variables may simultaneously belong to multiple groups, resulting in overlapping group structures. For example, in genomics, the same gene may participate in several
14

---

G1
G2
Go
f2:
Fpp Ri0e F11
Groups with Overlapping
Gi ias Fa
G2 Jaa Ha
G2 Zu
G4 Zts
C Sg.
Nonoverlapping Decomposition
Figure 2: Illustration of overlapping and disjoint group structures. The upper panel shows three pathways G1, G2, G3 with shared genes, leading to overlapping groups. The lower panel demonstrates the decomposition of these overlapping pathways into non-overlapping subgroups (G'...., G'²), which allows group-level dependencies to defined a consistent and interpretable manner.
biological pathways, different pathways share common genes. Such overlaps complicate the definition of group-level dependencies, as conventional methods typically assume disjoint group memberships.
To illustrate, consider three pathways denoted by G1. G2. and G3 as shown in Figure 2 Each contains five genes: G1 2, 3, 4, 5), G2 {4, 5, 6, 7, 7, G3 {7, 8,9, 10, 11). Here, G1 and G2 overlap on genes {). G2 and G3 overlap on {7,8), while G₁ and G; have no common elements. To resolve overlaps, decompose the original groups into disjoint subgroups shown in
Figure 2 G = {1, 23 G.2 (4,5). G. {6), (7,8), (9, 10. Accordinggn the original pathways can represented G G2 G2 G G3, G4), G3 = (G', G We denote l/ 0 whenever two disjoint subgroups Gj and Gy are conditionally independent. Based this decomposition, define the conditional dependency rules among the original pathways:
Non-overlapping groups: For instance, G1 and G3 have no overlap. If 0' 01 8/ = 82 0, then G1 and G3 are conditionally independent, as all variables in G1 are independent of those in G3. Otherwise, they are conditionally dependent.
Overlapping groups: For adjacent groups such G₁ and G2, if 6/ * or 0 they are conditionally dependent, since at least one unique variable in G1 depends on unique variable in G2. If instead 013 - 01 O, but least one of 0z 02 022 is nonzero, dependency arises through the common subgroup. Otherwise, G₁ and G2 are conditionally independent.
This decomposition preserves interpretability by mapping subgroup levell dependencies back to the original overlapping groups, while avoiding inflated false discoveries due to redundant memberships.
A.3 Proof of Statistical Properties
A.3.1 Proof of Theorem 3.2
Since h())) E F., there exists T(()) E H 1]s, ; Bo) such that /(x)) T(()). Refer to Corollary 3.1 in [5), there exists function (()) implemented by ReLU network with width W 38(
1 1)23s 1 - N [log2 (8N) and depth D 21 [(3) 1)2 M [8M)) + 2s such that
(k()) f()) < 198, (L8) (8)+ (NNM) 28/s 0 11
(13)
15

---

Knowing that (BOrT)² T
be rearranged as (u=li)xs 0 0 ((pp)-))))
have
(OoT")
uē. Thus, for any low-dimensional Relu network (() there always exists
high-dimensional sparse ReLU network f(x) satisfying f(ø) fax). This completes the proof.
A.3.2 Supporting Lemmas
ay 1. There exists Co 0 such that Da0y,th Cr. 2 co and |
Lemma A.I.
2 for all
k E Aj and a E H'
2. Denote (bxx) the vector obtained from a by setting [l:e = for all (i,k) E Ui) a,
Bj and o = for all (1,1") U (1) Pt. For a E H", (((x) also belongs to H.
Proof. By Theorem 3.2 we establish the uniform convergence:
/15(0-5)) faz (x (x=;))
as -
(14)
This implies that faz (z- + B[a, |x- serves as the unique minimizer of (a,). Consequently, the composite risk functional R(a") Σラ (af) attains its global minimum through coordinate- wise optimization. For any competing parameter arp E H., this construction ensures:
Ry())o;) Raja)
Vj l.....p.
(15)
The fundamental inequality
Ra(a)) min [xx - g(x-,))² S
min Ry(a,) RR ayeW,
(16)
holds with equality if and only fa, fa; almost surely, where g can any measurable function from the input space to the real numbers. Therefore, the identifiability condition
ag E H = fas; = laz: Vj =j=1... .....PP.
(17)
follows necessarily.
D
Assuming that such co exists, there exist oig E H- and k A, such that (k)) O. This means faos Saz does not depend on the related variable zs.
(OOCC ④ DrT aad
which is a contradiction.
2.
Denote Φ(a;) the sub-vector of ø(a)) corresponding to aj. Since a E H", we have Ja (ms, 1B5) Saz (MB,, a8;) Sa, (0, x5;)) Sola;) (mB,, 2gy), which implies ((a)
H.
Lemma A.2. Lerp o(log n), SD log S = O(+ then there exists e 0 such that
log n
(a) - -[(a) S
Va E W,
(18)
n
holds with probability at least 1 - 81 with 8, 4n (nt +1) logn
(32ent)
(logn)nd
e- 12
Proof. Let (a:;) (xij Sa, (mr,-)))²² and (x(x)) Σj (2:). For E Xm, denote
Rn(fz) 1 Σm (x;) and R(f) = (x))
Define
Q= = (a E xm 3Bo, ... Sa, E F s.t. [R(J) Rn)(aa) 28].
(19)
and
R R= (T. (r E xm x²; 3f. fa, E Fst. RR. R.((fa) 2 8)
2
(20)
16

---

Since {[R(f) - (f- and [R(f) Rx((fo) < 2 C ((.. - Rn(fo)) 22), we have
P(R) 2 P 3fa
".aa..Ja .....aaa Sap stt. [R()) - RR 2 2 and R(f) R, (n,( <2) E)
1 P fs 3fa. ... fa, E F, R(() R., (x.() a and Rn (n.).)l < </2) dP(r).
(21)
Noting that fa, maps into 0. 1] and Tij E [0, 1], we have (a;) E [0, 1] and l5(xx)) E [0.p]- Hoeffding's s inequality for bounded random variables shows that
P ((R(() (1)l S 2 exepp
nE²
1
exp
(22)
4p
2
for n > -.- Thus, P(Q) S 2R) for 2.
In order to bound IP)),, the technique of permutation and reduction to a finite class is useful. Denote that switch elements and n + i, for i in some subset of {1, 2, 2nj. A permutation is chosen uniformly at random from Tm- By Lemma 4.5 of M..
Tm the of permutations on 21, .... 2n} ....,
PRR) (ow R) S
max P(ox R). xEX2n
(23)
where the expectation is over æ and the probability is over permutations a.
Denote
{1)()) ((())=E; Eド (I, - ()= fa, E FF = 1.....pp N (E.F, be the covering number of F under the empirical (P()))) metric with ra- F) be
LF
Let
dius E. For a given sequence (Z1....,,Za) .In) ) E Xr" let Flz the subset T..", Define the uniform covering number
No (FFF. = max (l, L'. : : x
(24)
Suppose that x E 2 and let a minimal 5-cover for Crlz with respect to the L' metric. Pick G C F such that T Lglz and [ IT. Let X and E 2 that x = (r,s), 8). Suppose that f E F satisfies RR, - Rn(f))l 2 2 and there exists g E 9 such that D (x:)) (x,)) g. Then,
1 n
2n
RR, [g+))
(z)))
Σ 1x(x;)
71
72 imn+1
2n
Σ (z((x;) - (x(xx))) -
Σ
(1,(1;) - (1y(x,)) + +R,(ft) - Ra(fe)
1=1
i=
|R (f+) Rn, (fa)l
(z(() - (x,())))
1
Σ (lg(,,) - (x(())
i=--1
2n
2 [Rn(fr) Rn Rn,(fa)
(x,2x,) - (x(x)))
8/4..
(25)
17

---

Thus,
Plow
R) S 3g g:
1
Σ (l((0o(1)) (g(
2 8/4
a
S |9 I gEg
(10 (2)(1)) - (xFo(()))))
2
(26)
n
|g| max P gEG
Σ
(1(xz)) 1o (2n+tel 20/4
=²n
-
| exp
32p
where each & is independently and uniformly drawn from {-1,1}. Noting that
1
1 1pll))(P))
- (2:))
2n
1=1
1
Σ
Σ [1ig Say (zi [x0 Sa, (2i (113)
2n
(27)
1 2n
Σ (fa, (a,) + +,a - 2x15))
2n
1=1
1
Σ Iaa
fayl
can translate cover of Cflz into cover of the function space F, that is,
| = N (g. Lol,, (P.)) s N2n (G. s N2n 16p F.
(28)
For class of functions F, the pseudo dimension, denoted by (F) is natural measure of its complexity. According to [I. (F) is the largest integer m. for which there exists (.i,-.., Tm, V1,.-.,,.k Vm)) E xm x m such that for any (b1,..., ,,bm) {0, 1) there exists f E F such that Vi : f(x;) > = by = 1. Using Theorem 18.4 in [. we can give an upper bound on the covering number by Pd) Suppose that Pd)m(F) d, we have
2e d
Nº (E.FF, e(dd + 1)
(29)
Moreover, based on Theorems 3 and 6 in 12, the pseudo dimension of ReLU feedforward neural network space is bounded as
c c-SDlog (SDD) < PJim(F) < CSD log)
(30)
Thus,
P
22n
IP(Q)= < 4N2m
F. (Pn)
exp
16p
32p
(31)
32pe
DSDlog(S)
82n
$ _e (SD 1)P
exp
32p
Let p
o(log n),
O Jogn
and SD log S
O
it is obvious
that n
4 and
4e (SD log(S) + 1) (32pe)
PSD log(S)
exp
32p
2
4r (1 +1
log n (2ennt]
(log n)n
Ln log
0. This completes the proof.
18

---

Lemma A.3. There exists constant Ma, such that with probability at least 1 82, the empirical risk function (,) Ma, Lp-Liccchitz
Proof. For the empirical risk:
2
RR, (a) - (B)1 =
(xis - Sa, (x),)))) [zij - f3, (xi,-
J=1
-oo--faf (fa, (31, - fo, (T), (far (ot,))) + foz (x)1-s) - 2015)
fay (z) - f0, (31,-5))
/=1
fay (I1,-)) faz (Zx.-)) + [ (2o,-)) - faz (Ee +2|Ess)
Bl 3y11 . (2+ 2/s;y|)
/=1
i=1
2 n 2
|ag - -3y1
|e:g
n
(32)
For a zero-mean sub-Gaussian random variable Exj, we know that there exists a constant C₁ such that Eleyy s Cia. Applying Markov's inequality:
Σ E Bleisl
Eegl
Cba
S
(33)
n
Co
Choosing Co2 = =C20/82, we obtain
(a)(a) (B) S 20t) - -B1 BMM vplla B1,,
(34)
with probability at least 82-
A.3.3 Proof of Theorem 3.6
Define B, argmin | all. Let L(a) Ez= EE [
Eier 16zz and A
max(A, (A1, 22). Since L(xx) is Lipschitz, that is,
p-1
P
L(B,) L(ār)) Σ
+ Σ
10
Σ
Σ
10
ans
j=1
j=1 )=1 k=
LeF
8 p-1
| [s.i]]
Pnj
r(:Ak]
Cnny
+ Σ
18
Pam!!
(35)
j=1
I[ 1)
< Volp 1) [r, - -ra..l
100, - -0a.ll
2
C V p(p - 1) 1) - āall-
Combined with Lemma
A.2 we have
C2 |B âm = C2n H¹)² < R(ān) - R(B,)
logn
< 20
+ +X(L(B,) - L(aw))
(36)
log n
< 20
+ XC ppcppp - 1)113, - -á,l1
nt
19

---

holds with probability at least - 61. Applying Young's inequality,
v/(u-1)
AC - 1),1,, - -nnll 1 (c2u)1/v
11B (ann)nl) + v- / 20 100
plp 1)
A
(e2u)
02 1| - -nll + (AvP(p-1) (1ppp-1) w/(u-l)
$ 2
(37)
|| à C' ((w/plo-1) 1)
yielding
log
Let p ollog n), SD (S)
nt
O(n O(n Since v > 2, 1 < /
w-1
<
2, there exists C3 > such that d(órn; H²) S
1o
c3
Let K denote the inactive regularization components of L. Since is Lipschitz and K(ő(a)) = K(a), we have
P
A
Σ
A2
(n((((,)) R., R.,( M MK(()(B,)) K(án))
j=1
KEB
logn n
< 201
+ R()(B,)) R((am)
n
X[K K K(ā.)]
log n
< 2c1
XC p(p 1)113,, -á
ni
(38)
Similarly, since v > 2,
1
< 1, there exists
0 such that Σi=1 EKEB,
El Ever, 10,ll S c4 log
log n
This completes the proof.
A.3.4 Proof of Theorem 3.7
By Theorem and Lemma A.3 we have that for all (i, E Ui=1 x with probability 1 1- 81.
r (x)) and &c,,C, are bounded away from zero as + oc. Our analysis considers two connection structures: variable-level connection and group-level connection. The fundamental dependency principle requires that variables from statistically independent groups must exhibit no conditional
dependence. That is, 0c,,c. = 0 implics (+lee,(b]]] = Let M(a) = E P0
Eier
10vl and ( = max(Ci. (2). Thus, 16,17
p M(a*) Σ B
Σ
1
10
< 00,
(39)
r
Pep, 10
1=1
and
log
log
c2d(ān, H') < 201 ni + +((N(a") M((a)) S 2en ni + +(((a)).
(40)
Let Letç= 0 (n a there exists Gg > 0 such that d(ān, H²) c tr) 0 with
probability - 81. Thus, by Lemma A.I oc,.c. and (k)) are bounded away from zero for all k E Ay and large enough
To prove true negativity, we can separately prove that 1 (110,))) for all (jA)) E U3 x Bj and 81,r = 0 for all (1,4") e Ul_10 (4 x Pt. We establish the result by contradiction. Suppose there
20

---

exist some j and & E B, with 1 (10) ((()(()) component O, then we have (ā,( C2 p plog000 S R.,g. (g.). By Lemma there
* O. Define In the vector obtained from on by setting the ploy0001
exists Má, such that
F (
Rn((nn) & Md, vFlggn - -à,ll = MM vP / 1 C2 1hoooll S (g,
(41)
82- Since (1)
with probability at least
* O, we deduce that C2 p
M&a/P. This contradicts Theorem 3.6 which proves that for n large enough Ç2 1
S2 (log nogg))= - 2 2M vp. with probability at least 1 - 61. Thus, by Bonferroni inequal- ity, we have [1),, 0 for all & E B, with probability at least 1 - 81 82. Similarly, it can be inferred that ow for all l' Pi with probability least 51 62. This completes the proof.
A.4 Computation
A.4.1 Algorithm
Algorithm 1: Two-Stage Proximal Gradient Descent for NNBLNet
Input:
Data {a; learning rate y,, regularization parameters (i λI S2 A2, power ?, number of epochs T, tolerance
Output: Estimated parameters {0, T, (A 121)
Stage 1: Initial Estimation (non-adaptive)
Initialize 0, y. and A for all j and for i = to do
for j = 1 to p do
Compute predictions: Big = Say (s)
Compute gradients: Vr,, Ve. V 2 for all
Gradient step:
imp
-r,
@mmp 0 V
4 A v an
for all
Proximal update:
ni M
r
tmpp
|
Ow + (ign( max
7S2 0)
Store estimates r, i and compute adaptive weights: jijk
(2)
|r] 174e" WHy
+ 10
Stage 2: Adaptive Estimation
Reinitialize r, and 0, keep A from Stage (or reinitialize optionally) for i to 7 do for j= 1 to p do
Compute predictions and gradients in Stage
Gradient step and update A as before
Proximal update with adaptive weights:
tmmp
Ow & sign( max
( mmp
T/(2wl? (2)
21

---

A.4.2 Selection of the Regularization Parameters
Based on the convergence requirements specified in Theorems and 3.7 we established the parame- ter configuration A <* = en 1/8
1/8 for k 1, Computationally, n 17s yields values between 0.46 and 0.26 for sample sizes ranging from 500 to 50,000. To determine the optimal in Table we eval- uated Fi-score performance across values {0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 1, 1.5) at different sample sizes. This theoretically derived value achieved stable performance across varying n at c = 0.35. We therefore use c = 0.35 in practical implementations to maintain computational efficiency without compromising accuracy.
Table 5: Fl score of NNBLNet with different sample sizes and values of under linear and nonlinear scenarios
Sample size n
Sccnario
500
1000
2000
5000
10000
20000
0.05
0.659 (0.033)
0.666 (0.026)
0.673 (0.023)
0.674 (0.022)
0.670 (0.019)
0.676 (0.017)
0.10
0.695 (0.032)
0.704 (0.028)
0.713 (0.022)
(0.018)
0.719 (0.016)
0.717 (0.014)
0.15
0.720 (0.030)
0.729 (0.025)
0.731 (0.019)
0.738 (0.017)
0.735 (0.014)
0.742 (0.013)
0.20
0.738 (0.024)
0.749 (0.021)
(0.015)
0.760 (0.012)
0.758 (0.011)
0.25
0.756 (0.022)
0.763 (0.017)
0.769 (0.015)
(0.013)
0.771 (0.011)
0.773 (0.010)
Nonlinear
0.30
0.764 (0.017)
0.767 (0.015)
0.772 (0.013)
0.777 (0.011)
0.778 (0.010)
0.780 (0.009)
0.35
0.769 (0.019)
0.772 (0.014)
0.775 (0.012)
0.778 (0.009)
0.776 (0.009)
0.778 (0.008)
0.40
0.765 (0.020)
0.768 (0.017)
0.775 (0.013)
0.777 (0.011)
0.781 (0.008)
0.779 (0.008)
0.45
0.762 (0.021)
0.763 (0.019)
0.767 (0.014)
00.703)
0.773 (0.011)
0.771 (0.010)
0.50
0.748 (0.022)
0.755 (0.018)
0.757 (0.015)
0.762 (0.014)
0.761 (0.013)
0.765 (0.011)
1.00
0.708 (0.027)
0.716 (0.023)
0.724 (0.019)
0.729 (0.018)
0.728 (0.016)
0.731 (0.014)
1.50
0.670 (0.036)
0.677 (0.029)
0.688 (0.024)
0.695 (0.022)
0.693 (0.018)
0.698 (0.016)
0.05
0.645 (0.035)
0.652 (0.028)
0.660 (0.025)
0.661 (0.023)
0.657 (0.020)
0.664 (0.018)
0.10
0.682 (0.033)
0.691 (0.029)
0.701 (0.023)
0.704 (0.019)
0.707 (0.017)
0.705 (0.015)
0.15
0.707 (0.031)
0.716 (0.026)
0.718 (0.021)
0.726 (0.018)
0.723 (0.016)
0.730 (0.014)
0.20
0.724 (0.026)
0.735 (0.022)
0.741 (0.018)
0.747 (0.016)
0.747 (0.014)
0.746 (0.012)
0.25
0.742 (0.024)
0.751 (0.019)
0.757 (0.016)
0.761 (0.014)
0.760 (0.012)
0.763 (0.011)
Linear
0.30
0.751 (0.019)
0.755 (0.016)
0.760 (0.014)
0.765 (0.012)
0.766 (0.011)
0.768 (0.010)
0.35
0.760 (0.021)
0.765 (0.016)
0.769 (0.014)
0.772
0.771 (0.010)
0.773 (0.009)
0,40
0.752 (0.021)
0.755 (0.018)
(0.05)
(0.012)
0.767 (0.010)
0.766 (0.009)
0,45
0.746 (0.022)
0.749 (0.020)
(0.055)
0.757 (0.014)
0.760 (0.012)
0.759 (0.011)
0.50
0.732 (0.023)
0.741 (0.019)
(0.46)
0.749 (0.015)
0.748 (0.014)
0.752 (0.012)
1.00
0.695 (0.028)
0.703 (0.024)
0.711 (0.020)
(0.199
0.715 (0.017)
0.718 (0.015)
1.50
0.657 (0.037)
0.664 (0.030)
0.675 (0.025)
0.682 (0.023)
0.681 (0.019)
0.686 (0.017)
For applications requiring finer calibration, the regularization parameters may alternatively be selected via cross-validation. Specifically, avoid extensive grid search, we unified the four sparsification parameters as a single Ag based on theoretical analysis indicating their identical asymptotic order. Certainly, should that be necessary, we could alternatively assume distinct pa- rameters for the four components and conduct a grid search. We then implemented a five-fold cross-validation procedure: datasets were partitioned into training and validation sets, where mod- els were trained using objective functions (3) and 4 to obtain fa, (Zx.-)) for j = 1, 1,...,, ,p. The
validation loss L = = Σ (Ii.j Faa, (2i.-1)) 2 was evaluated across candidate An values (0.1. 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, 1.0), with the minimizer selected as optimal. Performance results are summarized in Table 6 This approach yielded slightly improved Fl scores compared to our original results but incurred higher computational complexity.
A.4.3 Selection of the Hyperparameters
Regarding the hyperparameters of the neural network, we evaluated different configurations using two simulated datasets. Specifically, we investigated the effect of the number of hidden layers: (2, 3, 5. 8) and the number of units per hidden layer: (25, 50, 100). The Fl-score results (summarized in Table 7 indicate minimal performance differences across configurations. To optimize computational efficiency while maintaining competitive performance, and in line with established practices in the field, we ultimately set the network architecture to 3 hidden layers with 50 units per layer.
22

---

Connections were introduced between each pair of adjacent groups, i.c., group 1 is connected to group 2, group 3 is connected to group 4, .... and group 9 is connected to group 10. For each such pair, nove-level links were established such that the i-th node in the first group was connected the Table 6: Performance of NNBLNet with the optimal λ identified by five-fold cross-validation across six datasets (mean values and standard deviations for synthetic datasets).
Dataset
Recall
Precision
F1
Nonlinear
0.878 (0.016)
0.702 (0.021)
0.776 (0.014)
Linear
0.893 (0.017)
0.695 (0.019)
0.781 (0.013)
Friendship
0.888
0.820
0.853
Co-authorship
0.732
0.698
0.715
BRCA
0.779
0.634
0.700
LUAD
0.658
0.545
0.596
Table 7: Fl score of NNBLNet with different hyperparameters under linear and nonlinear scenarios
Scenario Type
Layers
Units
Nonlinear
Linear
25
0.743(0.019)
0.739(0.018)
50
0.755(0.018)
0.751(0.017)
100
0.748(0.018)
0.742(0.018)
25
0.765(0.016)
0.760(0.015)
50
0.772(0.014)
0.769(0.013)
100
0.768(0.016)
0.765(0.015)
25
0.769(0.018)
0.765(0.017)
50
0.773(0.015)
0.770(0.014)
100
0.769(0.017)
0.766(0.017)
25
0.768(0.019)
0.764(0.018)
50
0.771(0.016)
0.768(0.016)
100
0.765(0.018)
0.762(0.018)
A.4.4 Computation Cost
All experiments were conducted on a workstation equipped with an Intel Core i7-800H Processor, an Nvidia Tesla A40 GPU, and 64GB of RAM. Table X compares the computational time and peak memory usage of our method against the baselines across varying values of n and To ensure fair comparison, all methods were executed the most computationally efficient method by significant margin. In contrast, NNBLNet, BGSL, and DeepGRNCS involve substantially computational costs, which escalate with the number of nodes. This higher cost is attributed to the iterative training of neural networks and, for BGSL, the additional overhead from MCMC sampling. Although NNBLNet slower than the anear-based methods, is faster than the other neural network-basse approach, DeepGRNCS. Overall, NNBLNet demonstrates favorable scalability.
single CPU core. The results indicate that Fair Glasso is
A.5 Experiment Details
A.5.1 Synthetic Data Setting
We simulated a dataset with n = 1000 samples and 100 variables. The 100 nodes (variables) were partitioned into L 10 groups (blocks), each containing 10 nodes. For each block, we generated a Aarabási-Albert network structure with maximum node degree of 4.
i-th node in the second group (e.g., I to 1 32 to 12 etc.).
Within each group, the features were generated using recursive formula inspired by the Barabási- Albert network property. First, nodes were sorted by degree in descending order: f1,..., For the node ji with the largest degree, Jijh was generated from standard normal distribution. For
23

---

Table 8: Average computation time (in minutes, the first value in parentheses) and peak memory usage (MB, the second value in parentheses) of different methods across various sample sizes (n) and dimensions (p)
p
NNBLNet
BGSL
Fair Glasso
DeepGRNCS
50
(9.5, 61.8)
(6.3, 61.8)
(0.0, 30.3)
(15.7, 240.4)
500
100
(21.8, 307.1)
(11.0, 90.3)
(0.0, 45.1)
(49.6, 1120.7)
200
(46.2, 709.5)
(14.5, 149.2)
(0.1.69.1) 69.1)
(100.9, 4490.8)
50
(12.2, 158.2)
(8.1, 65.7)
30.4,
(20.1, 338.8)
1000
100
(28.0, 308.9)
(14.1, 91.1)
(0.1,45.7)
(63.6, 1148.2)
200
(59.2, 719.6)
(18.6, 152.3)
(0.1, 70.2)
(129.4, 4602.9)
50
(15.6, 179.7)
(10.4, 70.1)
(0.2,
(25.7, 399.2)
2000
100
(35.8, 428.4)
(18.0, 98.0)
(0.3,
(79.0, 1652.6)
200
(75.8, 929.2)
(23.8, 171.9)
(0.4, 89.9)
(150.0, 6397.4)
50
(19.5, 215.8)
(13.0, 74.8)
(0.3,
(32.2, 528.7)
5000
100
(44.8, 498.9)
(22.6, 109.2)
(0.3,
(95.0, 1952.7)
200
(94.7, 1102.5)
(29.8, 182.5)
(0.7, 101.8)
(185.0, 7598.6)
50
(22.6, 241.0)
(15.0, 80.9)
(0.5, 38.5)
(37.2, 599.5)
10000
100
(51.8, 579.2)
(26.1, 114.7)
(0.8, 60.1)
(110.0, 2248.9)
200
(109.5, 1247.9)
(34.4, 193.3)
(1.1, 109.7)
(215.0, 8799.1)
50
(26.3, 257.6)
(16.8, 85.3)
(0.7, 39.8)
(42.4, 633.2)
20000
100
(57.3, 609.8)
(30.8, 127.2)
(1.0, 64.4)
(121.2, 2540.5)
200
(123.9, 1335.2)
(39.7, 208.8)
(1.4, 117.9)
(242.2, 9008.7)
each subsequent node (1 (1= 2,..., 10) with neighbors Nj. in the same group, the feature Zzj) was simulated as:
Ziji
E fx( Eggo
(42)
JeEN A<1
where for the nonlinear case: fx(() incorporated polynomial term 2², interaction term 0.5 0.2 ex and sin xr and for the linear case: fa (x) x with z ~ (2,, 1). The error term Eij, was generated from (00 0.01).
For each connected group pair (e.g., group and group 2), we first generated the features of the first group (e.g., group 1) using the procedure described above. Then, for the second group (e.g., group 2). we generated each node's feature using a bi-level formulation. Specifically, we first generated group-internal signal T'j EjeEN, fu (wxje). then incorporated signal from the connected
node in the previously generated group, g(x), along with group-connection coefficient 6* as follows:
~N(5, 1)
Zíj. lr (zly, + g(zaj.) Etji
(43)
where j' denotes the index of the connected node in the first group, and g(z) I3 for the nonlinear case and g(x) x for the linear case.
In the synthetic data generated by the aforementioned procedure, inter- group connections accounted for 37.0% the total network edges. To better visualize the connectivity patterns, we quantified the network's 's structural properties using two metrics. The intra- antra-group density of group C₃. is defined as the proportion of observed edges among all possible edges within the group:
IntraDensity(Ca) Eiej,eeE, Aij
[CW[[[a|-1)
Similarly, the inter- inter-group density between two distinet groups C₁ and Ct (k o 1) is given by the proportion of observed edges across the groups:
InterDensity(C,, Cj) Eiecr.sec Aij
1cal-105

---

For the synthetic network, the mean intra- and inter-group densities were 37.8% and 10.0%, re- spectively. These metrics help characterize the strength of connections within and between groups. providing quantitative basis for analyzing the network structure.
For both linear and nonlinear patterns, the performance was assessed over 100 independent simulation replicates, with results for recall, precision, and Fl-score summarized by their mean and standard deviation.
A.5.2 Real-world Datasets
Table summarizes the key features of the real-world datasets. Additional details are provided below.
Table 9: Summary of real-world datasets.
Dataset
Nodes (No.)
Edges (No.)
Samples (No.)
Groups (No.)
Friendship
Students (311)
1009
Interactions (47127)
Gender (2)
Co-authorship
Authors (130)
525
Keywords (1903)
Publication type (6)
BRCA
mRNA (73)
763
Patients (1099)
Pathway (4)
LUAD
CNA (98)
700
Patients (507)
Pathway (13)
Friendship. This network dataset captures social interactions among students across nine classes a high school in Marseille, France, recorded over five consecutive days in December 2013, Following standard practice in contact network analysis, we constructed a ground toound-truth graph where nodes correspond to students and weighted edges reflect aggregated interaction frequencies. In line with [25, node attributes were assigned based gender, and the signals were generated by grouping the terns
interactions into sets of four. The Friendship dataset is available at www.aaddfiiiiindddeeeeeeeeeeeeee
Co-authorship. This network dataset originates from ACM conference proceedings and includes 17,431 unique authors, 122,499 publications, and 1,903 technical keywords. We focused on representative subset of authors, where nodes correspond to individual researchers. Demographic attributes were assigned according to authors' predominant conference categories, determined by their maximum publication frequency. The ground-truth network was constructed through co- authorship detection: edges were added between authors who co-published least paper. To generate network signals, we quantified authors' keyword usage patterns by calculating normalized frequencies of specific technical terms across their publications. The Co-authorship dataset is available at atttttss: /tttpp::/dd.cee acm. org/
BRCA. The BRCA dataset is derived from The Cancer Genome Atlas (TCGA) and comprises mRNA gene expression profiles from 1,099 breast cancer patients. We selected 73 genes involved in four key biological pathways: B Cell Receptor Complexes, Caspase Cascade, G1 And S Phases, and MMP Cytokine Connection, based on prior domain knowledge and pathway annotations from the KEGG database [T6.. Nodes in the network represent these genes, and group structure is defined by their pathway annotations [42]. A biologically grounded reference network was constructed using curated interaction data from the STRING database (35), which integrates multiple evidence sources such as experimental data, co-expression, and pathway information. An undirected edge was placed between two genes if high- nigh-confidence interaction reported in STRING, reflecting known regulatory or functional associations. For network estimation, each patient's 's expression profile was treated an input signal across the 73 genes.
LUAD. The LUAD dataset is also sourced from TCGA and contains copy number alteration (CNA) profiles for 507 lung adenocarcinoma patients. We focused on 98 CNAs implicated in 13 distinct biological pathways (see Figure 3 for details), with pathway annotations obtained from the KEGG database, offering more complex grouping structure compared to BRCA. Each CNA is treated as node, and group membership is determined by pathway assignment. As with BRCA, the reference network was constructed using functional interaction information from the STRING database. Each patient's CNA profile was used as input for structure learning, with the goal of recovering sparse and modular dependencies.
BRCA and LUAD expression data were obtained from the R package edgsr, pathway information from the KEGG database was obtained using msigdbr, and interaction information from STRING was
25

---

obtained via STRINGdb. The structures of the four benchmark networks reveal distinct connectivity patterns. Specifically, the proportion of inter-group connections is 47.9% (Friendship), 36.2% (Co- authorship), 60.8% (BRCA), and 74.4% (LUADD, while the corresponding intra- and inter-group density pairs are (37.8%, 10.0%); (2.5%, 2.0%); (11.1% 1,, 3.2%); and (53.7%, 20.9%). These metrics collectively highlight the structural heterogeneity across the datasets.
A.6 Sensitive Analysis
A.6.1
Generalization Capability
To evaluate the generalization capability of the proposed method, assessed the stability of its network estimates. We repeatedly drew 90% subsets of the data from two synthetic and four real- world datasets and re-estimated the networks 100 times. The Jaccard index between each re-estimated network and the original full-dataset network was then calculated. The resulting indices (mean and standard deviation) were as follows: Nonlinear: 0.913 (0.052), Linear: 0.932 (0.047), Friendship: 0.889 (0.061), fo-authorship: 0.921 (0.055), BRAC: 0.902 (0.058), and LUAD: 0.874 (0.064). These consistently high values indicate that the proposed method yields stable network structures across diverse datasets, confirming its strong generalization capability.
A.6.2 Performance across Varying Sample Sizes and Group Sizes
To assess the robustness of method sample and group size specifications, we evaluated its performance under varying sample sizes (n = 500, 1000, 2000, 5000) with a fixed group size of ten, and under varying group sizes (ranging from 5-10, 5-20, and 5-30 individuals) with fixed sample size of 1000. For groups containing 20 or 30 members, the network structure was simplified to a star module with node as the hub connected setting (with a ratio of group size 1:1), these configurations reflected increasing levels of group size imbalance, with approximate ratios of 2:1, 4:1 and 6:1, respectively. The last case represents the most pronounced disparity, where the largest group is six times larger than the smallest.
all other nodes. Based on the balanced baseline
Results are summarized in Tables 10 and I Consistent with expectations, larger sample sizes enhanced the network recovery accuracy of all methods; nevertheless, NNBLNet consistently demon- strated superior performance. According Ta NNBLNet maintained highly stable performance across moderate group size imbalances, with only slight decrease in recall, precision, and F1 score observed under the most extreme 6:1 condition.
Table 10: Performance comparison of different methods for the nonlinear synthetic network with varying sample sizes: Mean (SD) over 100 replicates.
Method
Recall
Precision
Fl
NNBLNet
0.861 (0.022)
0.679 (0.027)
0.769 (0.019)
Fair Glasso
0.756 (0.021)
0.639 (0.026)
0.693 (0.020)
500
BGSL
0.764 (0.022)
0.631 (0.024)
0.692 (0.019)
DeepGR
0.782 (0.021)
0.662 (0.023)
0.717 (0.018)
NNBLNet
0.872 (0.016)
0.693 (0.022)
0.772 (0.014)
Fair Glasso
0.779 (0.014)
0.653 (0.021)
0.710 (0.017)
1000
BGSL
0.790 (0.017)
0.644 (0.018)
0.709 (0.013)
DeepGRNCS
0.809 (0.016)
0.681 (0.018)
0.731 (0.017)
NNBLNet
0.879 (0.014)
0.707 (0.019)
0.775 (0.012)
Fair Glasso
0.787 (0.013)
0.660
0.717 (0.015)
2000
BGSL
0.795 (0.015)
0.652 (0.017)
0.717 (0.014)
DeepGRNCS
0.822 (0.014)
0.695 (0.017)
0.752 (0.015)
NNBLNet
0.887 (0.011)
0.716 (0.015)
0.778 (0.009)
Fair Glasso
0.795 (0.011)
0.667 (0.015)
0.725 (0.012)
5000
BGSL.
0.802 (0.012)
0.660 (0.014)
0.724 (0.011)
DeepGRNCS
0.838 (0.011)
0.708 (0.014)
0.767 (0.012)
A.6.3 Performance when Group Label Misclassification
Even when genuine group labels (e-g., pathway information) are available, misclassification may still occur and impair model performance. To assess the impact of such label errors, we randomly
26

---

Table 11: Performance comparison of different methods for the nonlinear synthetic network with varying group imbalance ratios: Mean (SD) over 100 replicates.
Ratio
Method
Recall
Precision
F1
NNBLNet
0.872 (0.016)
0.693 (0.022)
0.772 (0.014)
Fair Glasso
0.779 (0.014)
0.653 (0.021)
0.710 (0.017)
1:1
BGSL.
0.790 (0.017)
0.644 (0.018)
0.709 (0.013)
DeepGRNCS
0.809(0.016)
0.731(0.017) (.6810)
NNBLNett
0.874 (0.017)
0.691 (0.023)
0.771 (0.015)
Fair Glasso
0.773 (0.015)
0.649 (0.022)
0.706 (0.018)
2:1
BGSL
0.782 (0.018)
0.640 (0.020)
0.704 (0.015)
DeepGRNCS
0.802 (0.017)
16600000019)
0.730 (0.016)
NNBLNet
0.869 (0.018)
0.688 (0.024)
0.769 (0.016)
Fair Glasso
0.766 (0.016)
0.642 (0.023)
0.701 (0.018)
4:1
BGSL
0.770 (0.019)
(0.021)
0.698 (0.016)
DeepGRNCS
0.812 (0.017)
0.684
0.734 (0.016)
NNBLNet
0.854 (0.020)
0.675 (0.026)
0.753 (0.017)
Fair Glasso
0.755 (0.018)
0.634 (0.024)
0.692 (0.019)
6:1
BGSL
0.756 (0.020)
0.627 (0.023)
0.687 (0.017)
DeepGRNCS
0.808 (0.018)
0.677 (0.020)
0.727 (0.017)
scrambled node group labels in simulated data under nonlinear setting, using misclassification rates of 10%, 20%, and 30%, with each scenario repeated 100 times. We evaluated three group-aware Fair Glasso, and BGSL-under these conditions. As summarized in Table [2] all three methods exhibit performance degradation label error increases. Nonetheless, NNBLNet maintains relatively higher F1 scores and consistently outperforms both BGSL and Fair Glasso, demonstrating its robustness to moderate levels of label misclassification in real-world applications.
methods-NNEE
Table 12: Performance comparison of different methods for the nonlinear synthetic network with varying levels of group label misclassification: Mean (SD) over 100 replicates.
Error rate
Method
Recall
Precision
F1
NeBLNet
0.872 (0.016)
0.693 (0.022)
0.772 (0.014)
0%
Fair Glasso
0.779 (0.014)
0.653 (0.021)
0.710 (0.017)
BGSL.
0.790 (0.017)
0.644 (0.018)
0.709 (0.013)
NeBLNet
0.821 (0.022)
0.598 (0.029)
0.691 (0.019)
10%
Fair Glasso
0.756 (0.023)
0.564 (0.031)
0.662 (0.020)
BGSL.
0.742 (0.024)
0.553 (0.032)
0.651 (0.021)
NNBLNet
0.765 (0.026)
0.541 (0.034)
0.633 (0.022)
20%
Fair Glasso
0.702 (0.027)
0.509 (0.035)
0.602 (0.023)
BGSL.
0.688 (0.028)
0.496 (0.036)
0.588 (0.024)
NNBLNet
0.702 (0.031)
0.489 (0.039)
0.574 (0.026)
30%
Fair Glasso
0.648 (0.032)
0.455 (0.040)
0.540 (0.027)
BGSL
0.632 (0.033)
0.442 (0.041)
0.526 (0.028)
A.7 Downstream Analysis of the LUAD dataset
Based on the low-level network inferred by NNBLNet, the top three hub genes with the highest degrees are TP53, JUN, and CD44. All three are strongly supported by existing literature in the context of lung adenocarcinoma (LUAD) and lung cancer biology. Specifically, TP53 is one of the most frequently mutated tumor suppressor genes in LUAD, with mutation rates often exceeding 40-50%. Its loss mutation contributes to genomic instability, aggressive tumor behavior, and poor prognosis in lung cancer 3 JUN has been identified as a gene associated with response to PD-1 blockade, suggesting its potential biomarker for immunotherapy efficacy in non-small cell lung cancer (37. CD44, a recognized cell surface glycoprotein and cancer stem cell marker, is commonly used to identify stem-like subpopulations in lung cancer [12]
Based on NNBLNet analysis, the most significant pathway-level connection was observed be- tween "MMP Cytokine Connection" and "IL4 Receptor in B Lymphocytes", with extensive gene
27

---

connectivity indicating functional cross-talk. IL-4 experimentally regulates MMP expression; for instance, it suppresses IL- L--I-indueed MMP-3 transcription in human fibroblasts by inhibiting AP- I promoter binding [31]. In LUAD,, MMPs are markedly overexpressed and contribute to ECM remodeling, tumor invasion, and metastasis 43,, while II.-4 signaling modulates immune activity and tumor-microente micreevimeetenniineeennttt crosstalk [7]. This supports a model wherein I1.4//IL-4R signaling influences ECM dynamics in LUAD by modulating MMP pathway activity, potentially via transcrip- tional regulation. As shown above, the analysis further validates that the NNBLNet method can yield biologically meaningful networks.
7 COLL RCCEPTOR COMPICEDO
F13 3ANILISS
VVP OYTOMINC CONNOCTION
TRVA ROCOPTOR
SIGVALING PAT NAY
GI AMD IVSCS
PROGRVAINTT COLL DCXTI
CASPASE CASCADE
C2 AND THASES
FTEI FATMWAY
CHOMOTANS
LYPHDOYTER
CANDADE CVGUP EXPR
Figure 3: Network reconstruction of the LUAD dataset using NNBLNet. Nodes are color-coded by annotated pathways, and the top three hub genes are labeled.
28

---

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

---

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

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: The supplementary material includes the code used for the experiments presented in this paper. Links to all publicly available datasets, as well as instructions for generating the synthetic datasets used, are provided in Appendix A.5
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.ca
public//uiids//ooees Seaerideee//odSbmiiiccc wwt.miissssooocoiooooo for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run reproduce the results. See the NeurlPS code and data submission guidelines (https: ComeeSumieiiicc / "nipp.cccpubblice / ssi onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones omitted from the script and why.
At submission time, preserve anonymity, the authors should release anonymized versions (if' applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The hyperparameter choices and optimization algorithm are described in Appendix A.4
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core the paper to a level detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [No]
Justification: The results of the synthetic simulations reported Sections 5.1 are averaged over 100 independent trials, with both the mean and standard error computed. However, we observed consistently low variance across these realizations, and thus omit error bars the plots for visual clarity. In Section 5.2]. the experiments are conducted on real-world datasets using all available samples and statistical significance tests are not applicable.
Guidelines:
The answer NA means that the paper does not include experiments.
31

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative eтror rates).
show in tables or
error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: We provide a summary the compute resources used for the experiments in Appendix A.4.4
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary failed experiments that didn't make into the paper).
9. Code ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics ht
iel ines?
Answer: [Yes]
Justification:
The paper conforms with the NeurlPS Code of Ethics in every respect.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g.. there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The broader impacts of this work are discussed in Section 6.2
32

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path any negative applications, the authors should point it out, For example, it is legitimate to point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses no such risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
The datasets used in the experiments are properly cited.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (c.g., CC-BY 4.0) should be included for each asset.
33

---

For scraped data from a particular source (e.g., website). the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, Vatasets.papersstt paperswi thcode. comdetts has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmodel as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
34

---

Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if' applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, or formatting purposes and does not impact the core methodology, scientific rigorousness, originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does not involve LLMS as any important, original, non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non non-standard components.
Please refer to our LLM policy for what should or should not be described.
35


Reviews:
[{"review_id": "30SNm1aZ0l", "reviewer": "Reviewer_cx7U", "summary": "This paper proposes DeepBLNet for recovering bilevel graphical structures from data.\nThe method is designed to learn both a high-level adjacency structure and a low-level structure simultaneously, from observed data.\nThe authors model this task via a bilevel architectural decomposition.\nThe paper provides a consistency theorem that claims the proposed architecture can recover the true bilevel structure asymptotically under certain assumptions.\nExperimental validation is primarily conducted on synthetic datasets with known structure.", "questions": "1. How sensitive is your model to group assignment errors or group size imbalance?\n2. What is the inductive bias that makes DeepBLNet better than training a flat structure learner on a large graph?\n3. Do you plan to apply this method in real application domains, and if so, which ones?", "limitations": "yes", "rating": 4, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Pros:\n- Novel Problem Formulation: The paper addresses a new variant of structure learning: recovering bilevel graph structures.\n- Methodological Design: The proposed DeepBLNet is a well-designed neural network architecture that handles group-level and intra-group adjacency jointly\n- Theoretical Guarantee: The authors establish a consistency theorem for structure recovery. While based on idealized assumptions, this provides a valuable theoretical anchor.\n- Modular Framework: The approach is likelihood-model agnostic (Gaussian, Ising, etc.), making it flexible across domains.\n\nCons:\n- The consistency result is solid in isolation, but it depends on strong assumptions: perfect optimizer, sufficient model expressiveness, asymptotic data. There's no discussion of finite-sample behavior, generalization, or network capacity.\n- While the paper includes some experiments on real-world datasets (e.g., neuroimaging and gene expression), the evaluations are mostly structural and unsupervised. There is limited benchmarking against other models or external validation of whether the recovered graphs are meaningful for real tasks.\n- The paper lacks an architectural ablation study. In particular, it would be useful to understand whether the bilevel separation (group-level and intra-group adjacency prediction) actually improves performance, compared to a flat or shared model. Some ablation or alternative design comparison would clarify the necessity of the proposed structure.", "was_revised": true, "final_justification": "The authors' additional clarifications on finite-sample analysis, ground-truth validation for real datasets, and planned ablations alleviate many of my earlier concerns. While I still believe that some evidence (e.g., finite-sample robustness, downstream biological validation) is missing in the current version and relies on promised additions, I appreciate the detailed responses and the clear plan to strengthen the paper. I will keep my weak accept recommendation, as the paper is novel and relevant. I expect the authors to follow through on their promised ablation experiments in the final version to strengthen the contribution.", "rebuttal": "Dear Reviewer,\n\nThank you for your thorough review and invaluable feedback on our manuscript. We fully concur with your critical suggestions regarding finite-sample analysis, validation on real-world datasets, and ablation studies. Below is our point-by-point response:\n\n1. Response to Comment 1:\n\nOur theoretical assumptions (sparsity, smoothness, boundedness, population risk) align with established DNN research [Chen et al. 2021; Jiao et al. 2023; Liang et al. 2018] and reflect practical realities: Gene networks exhibit localized interactions (sparsity), gradual nonlinear effects (smoothness), and physical constraints (boundedness). We will enhance intuitive explanations of these foundations in revised contexts.\n\nWe fully recognize the importance of finite-sample analysis. In the paper, we validated the model’s performance under limited samples through two simulation experiments and four real-world datasets. In the revised paper, we will add simulation experiments with various sample sizes and examine the network recovery capacity via the Recall, Precision, and F1 scores. Moreover, as an unsupervised model, we cannot validate generalization ability on test sets in the manner of supervised learning. Instead, we will attempt to characterize the model’s generalization capability through the stability of network estimation. Specifically, by repeatedly sampling 90\\% of the data subsets and re-estimating networks, we will compute the overlap ratio of identified edges against the original full-data network. Similar to existing DNN research and unlike traditional linear models, our method requires a certain sample size. In addition, with increasing sample size, DeepBLNet tends to have higher network recovery accuracy and stability.\n\n2. Response to Comment 2:\n\nIn the current work, we conducted two sets of simulation experiments involving both linear and nonlinear relationships between variables, along with four real-world data analyses encompassing social and biological network construction. In the simulated data, the true network structure is known based on the data-generating mechanisms. For the real-world datasets, well-established ground-truth networks exist. Specifically, for the Friendship and Co-authorship datasets, the former was constructed using student contact frequency while the latter was built based on co-publication records between authors. Both methodologies are widely adopted in existing literature [Navarro et al. (2024)]. The remaining two biological networks were derived from the STRING database, a universally recognized resource for characterizing gene-gene interactions [Szklarczyk et al. (2025)]. By comparing the edges in our reconstructed networks against these ground-truth references, we can validate the method's effectiveness. \n\nIn the revised version, we will provide more detailed explanations of these benchmark networks. Additionally, we will perform further downstream analyses on the estimated networks, including hub identification and community detection. Specifically, for the constructed biological networks, we will examine the biomedical significance of identified hub genes and analyze shared biological functions within detected communities, thereby informing clinical strategies for disease diagnosis and prognosis.\n\n3. Response to Comment 3:\n\nWe fully acknowledge the importance of ablation studies. In the original manuscript (Tables 1-3), we designed DeepNet as a baseline, which disregards the bi-level structure. Our proposed DeepBLNet consistently outperforms DeepNet in network reconstruction accuracy. In the revised version, we will restructure the paper to include a ``Ablation Study'' section. Building upon the existing comparison between DeepBLNet and baseline DeepNet, we will introduce a new comparator method that estimates group-wise networks using group labels only. This tripartite framework will be augmented with sensitivity experiments across varying proportions of inter-group edges (10\\%–70\\% of total edges), collectively demonstrating the distinct advantages conferred by the bi-level architecture beyond group information alone.\n\n4. Response to Question 1:\n\nWe acknowledge that group assignment errors and group size imbalance are unavoidable challenges in practical applications. DeepBLNet is expected to exhibit a certain degree of robustness in facilitating group assignment errors and group size imbalance. This is because the inter-group connection strength is estimated based on aggregated signals from all nodes within a group. The misassignment of a minority of nodes will be averaged by the correctly assigned majority within the group. Furthermore, the inter-group penalty term treats all groups equally, preventing smaller groups from being overlooked. Furthermore, it should be noted that there is a certain level of group size imbalance in the real-world datasets analyzed in the paper. Specifically, the (smallest group size, largest group size) for the Friendship, Co-authorship, BRCA, and LUAD datasets are (141, 170), (2, 46), (15, 24), (5, 28), respectively, demonstrating some imbalance. Our method has delivered satisfactory results in all these scenarios.\n\nIn the revised paper, we will introduce new sensitivity experiments to quantify the impact of group size imbalance and grouping errors. Specifically, we will randomly misassign the group labels of a portion of nodes in simulated data at group error rates of 10\\%, 20\\%, and 30\\%. We will also consider group size settings of (smallest group size, largest group size) = (5, 10), (5, 20), and (5, 30), corresponding to imbalance ratios of 2:1, 4:1, and 6:1, respectively. Preliminary exploration indicates that DeepBLNet exhibits a certain degree of robustness. \n\n5. Response to Question 2:\n\nThe core inductive bias of DeepBLNet posits that dependency relationships in complex systems exhibit a hierarchical structure. This bias is implemented through two key mechanisms: intra-group information sharing, where variables within the same group share latent functions (e.g., biological pathways), enabling the amplification of weak connection signals both among intra-group variables and across different groups through aggregation; and inter-group sparse transmission, whereby cross-group connections must be regulated via inter-group “switches” $(\\theta_{ll'})$, effectively preventing noise propagation between unrelated groups. Such hierarchical structuring finds strong empirical grounding in real-world complex systems. For instance, in biological systems, genes are grouped into pathways based on functional similarity, resulting in inherently layered dependency networks encompassing both higher-order pathway-to-pathway interactions and fundamental-level gene-to-gene relationships. The advantages of this bilevel network architecture are well-established in previous Gaussian graphical models designed for linear dependencies [Shan et al. (2020)]. Building upon this foundation, our work extends this conceptual framework into the nonlinear domain using deep neural networks (DNNs), thereby enabling more effective characterization of complex dependency relationships in real-world networks. In the revised manuscript, Section 5 will provide expanded treatment of this architectural approach.\n\n6. Response to Question 3:\n\nWe sincerely appreciate your concern about the practical application of our method. To substantiate its utility, we have systematically implemented DeepBLNet in four real-world scenarios: Friendship, Co-authorship, BRCA, and LUAD, spanning two critical domains: social networks (Friendship and Co-authorship) and biological networks (BRCA and LUAD). The constructed networks were rigorously validated against widely accepted ground-truth networks, confirming their accuracy. These applications not only demonstrate the method's versatility but also highlight the core value of our bilevel network modeling in these fields. While the current study focuses on these domains, the hierarchical modeling principle underlying DeepBLNet is broadly applicable to other complex systems involving group structure, such as financial networks, sensor networks, or organizational dynamics. In addition, the resulting networks readily enable downstream analyses, such as hub node identification and community detection. Specifically for biological networks, the identified hub nodes and functionally coherent gene communities within detected modules may provide critical insights into disease pathogenesis and progression mechanisms. Similar analytical approaches have been successfully employed in prior studies [Wu et al. (2023)]. In the revised manuscript, relevant downstream analyzes and discussions will be added.\n\n\nReference:\n\n[1] Chen, Y., et al. (2021). Nonlinear variable selection via deep neural networks. Journal of Computational and Graphical Statistics,30(2), 484-492.\n\n[2] Jiao, Y., et al. (2023). Deep nonparametric regression on approximate manifolds: Nonasymptotic error bounds with polynomial prefactors. The Annals of Statistics, 51(2), 691-716.\n\n[3] Liang, F., et al. (2018). Bayesian neural networks for selection of drug sensitive genes. Journal of the American Statistical Association, 113(523), 955-972.\n\n[4] Navarro, M., et al. (2024). Fair glasso: Estimating fair graphical models with unbiased statistical behavior. Advances in Neural Information Processing Systems, 37, 139589-139620.\n\n[5] Shan, L., et al. (2020). Joint estimation of the two-level gaussian graphical models across multiple classes. Journal of Computational and Graphical Statistics, 29(3), 562-579.\n\n[6] Szklarczyk, D., et al. (2025). The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research, 53(D1), D730-D737.\n\n[7] Wu, Y., et al. (2023). Identification of the hub genes in polycystic ovary syndrome based on disease‐associated molecule network. The FASEB Journal, 37(7), e23056."}, {"review_id": "4mLUmpnemI", "reviewer": "Reviewer_8Kg1", "summary": "This paper proposes a bi-level network estimation model that learns dependencies between variables at both the cluster level and the individual node level, connecting the two via a selection-consistent mechanism. The method outperforms baseline approaches on real-world and synthetic networks with nonlinear signals, and achieves comparable results on synthetic networks with linear signals.", "questions": "1. *Clarifying the Utility of the Bi-Level Architecture*. The paper uses DeepNet as a baseline to demonstrate the advantage of the bi-level architecture. However, if I understand correctly, DeepNet does not utilize the label information, which raises a question: to what extent does the observed improvement stem from the architecture itself, versus the use of additional label information?  This comparison will be more convincing if the authors willing to consider the following. Suppose the inter-cluster connections are sparse, thereby naively predicting them to be absent will only cause a small decrease in accuracy. Then, evaluate DeepNet where training is done within each cluster. In this case, DeepNet can leverage the label information to a partial extend. If DeepBLNet still outperforms this modified DeepNet, it would provide stronger evidence for the value of the bi-level design rather than just label availability.  Additionally, reporting the within-cluster and between-cluster edge densities (including synthetic networks) would help clarify how much signal is available at each level of structure, and how much the model is leveraging inter- vs. intra-cluster information.\n\n2. *Group level sparsity*. In Equation (3), the Lasso penalty on theta implicitly assumes a priori that the inter-group connections are sparse. However, this may not hold in all settings. In some cases, θ could be dense, but with lower (or higher) average weights compared to within-group connections—for instance, in the contextual stochastic block model. While sparsity is a reasonable assumption at the individual edge level, it is not necessarily justified at the group level. The paper would benefit from a clearer discussion of when and why this sparsity assumption is appropriate. This is also related to the last question of whether a prior knowledge on sparse inter-group connections contribute mostly to the performance improvement.  \n\n3. *Scale-free terminology*. Referring to a network with only 10 nodes and a maximum degree of 4 as \"scale-free\" is somewhat misleading. The concept of a scale-free network typically involves a wide range of scales (i.e., node degrees spanning multiple orders of magnitude). With such a small network, it is difficult to meaningfully assess whether the degree distribution follows a power-law or exhibits true scale-free behavior.", "limitations": "Yes", "rating": 4, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "**Strengths**\n\nThe idea of network estimation from a multiscale perspective is novel and compelling, especially considering that many real-world networks exhibit hierarchical or community structures. While multiscale modeling is well established in fields such as PDE-based dynamics, it is less explored in network estimation, making this approach interesting.\n\n**Weaknesses**\n\n1. As acknowledged by the authors, the proposed method relies on known node labels that reflect the community structure. This assumption could limit the method’s applicability. It would be helpful if the authors discussed typical real-world scenarios where such labels are available and justified their practicality. In addition, it is not entirely clear how much of the performance gain is attributable to the inter-class sparsity information provided by the labels, versus the bi-level architectural design of the proposed method. See the Questions section for further discussion.\n\n2. How the parameters $\\zeta$ and $\\gamma$ are chosen should be discussed, as they directly controls the level of sparsity, and can reflect strong prior knowledge of the network. For the network estimation problem, it seems there is no such a validation set in determine these parameters. Also it will be helpful if the dependency on these parameters can be represented.", "was_revised": true, "final_justification": "My main concern is that the proposed method relies heavily on prior knowledge and careful control of sparsity. However, the additional experiments presented during the rebuttal demonstrate that the algorithm maintains a certain level of stability under these conditions.", "rebuttal": "Dear Reviewer,\n\nWe sincerely thank you for your profound insights. We have systematically addressed core concerns regarding group label universality, sparsity parameter selection, and contribution decoupling of hierarchical structures. Below is our point-by-point response:\n\n1. Response to Comment 1:\n\nThe widespread availability of valid group labels in practical applications is well-established. For instance, in genomics, the majority of human genes possess annotated pathway labels indicating shared functional characteristics, accessible through biological databases such as KEGG. Similarly, public social datasets commonly include group division labels, as evidenced by gender classifications among students or research domain categorizations for academic authors. All four real-world datasets analyzed in Section 4.3 incorporate such group labels, which have been extensively utilized in prior research [Navarro et al. (2024), Szklarczyk et al. (2025), Lingjærde et al. (2021)]. We will expand the discussion regarding these group labels in Appendix A.5.2. Furthermore, even when prior labels are unavailable, unsupervised methods like spectral clustering can be employed to estimate group memberships, thereby ensuring methodological universality. Related considerations will be addressed in Section 5.\n\nRegarding the of contributions our bi-level architecture design, we will incorporate an ablation study subsection in the revised manuscript. The specific experimental design and outcomes will be detailed in the response to Q1.\n\n2. Response to Comment 2:\n\nGiven the unsupervised nature of our problem, which makes conventional train-validation splitting not applicable for parameter selection via cross-validation, we established the parameter configuration $\\lambda_k=\\eta_k = c n^{-1/8}$ for $k=1,2$ based on the convergence requirements specified in Theorems 3.6 and 3.7. Here, $c$ denotes a tunable constant and $n$ represents sample size. This formulation aligns with established practices in prior sparsity-inducing models [Liu et al. (2017)]. Computationally, $n^{-1/8}$ yields values between 0.46 and 0.26 for sample sizes ranging from 500 to 50,000. To determine the optimal $c$, we evaluated F1-score performance across values $\\{0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 1, 1.5\\}$ at different sample sizes. The results showed stable performance across a substantial range of $c$ values, with $c = 0.35$ positioned centrally within this stable region. We therefore use $c = 0.35$ in practical implementations to maintain computational efficiency without compromising on accuracy. Supporting results will be presented in the supplementary documentation of the revised manuscript.\n\n3. Response to Question 1:\n\nFollowing your suggestion, in the revised manuscript, we will incorporate a new ablation study subsection. This experiment employs a three-way comparative framework: alongside the original DeepBLNet and one baseline DeepNet, we introduce Modified DeepNet where groups are trained independently. Additionally, we conduct group structure sensitivity analyses across varying proportions of true inter-group connections (10\\%-70\\% of total edges) to isolate performance advantages attributable specifically to bi-level structure versus group label information. Preliminary exploration indicates DeepBLNet has superior performance, because even when inter-group connections are sparse, the relative proportion of inter-group edges remains significant given the equally sparse intra-group connectivity. Neglecting these inter-group connections would substantially compromise estimation accuracy. Crucially, real-world systems inherently require modeling inter-group connections: as evidenced by social contacts between gender groups in Friendship networks, interdisciplinary collaborations like statistician-biologist partnerships in Co-authorship data, and biologically established cross-pathway gene interactions in BRCA and LUAD datasets. Therefore, approaches that exclusively leverage intra-group connections while ignoring inter-group relationships constitute fundamentally inadequate modeling strategies.\n\nIn our simulation experiments and four real-world data analyses, the proportions of inter-group connections relative to total network edges are 37.0\\%  (synthetic network), 47.9\\%(Friendship), 36.2\\%(Co-authorship), 60.8\\%(BRCA), and 74.4\\%(LUAD). In addition, the mean intra-cluster densities and inter-cluster densities for synthetic network and Friendship, Co-authorship, BRCA, and LUAD networks are (37.8\\%,10.0\\%), (2.5\\%,2.0\\%), (11.1\\%,3.2\\%), (53.7\\%,20.9\\%), (44.0\\%,12.2\\%), respectively. The intra-cluster density of a group $C_k$ is defined as the proportion of observed edges among all possible edges within the group:\n$\n\\text{IntraDensity}(C_k) = \\frac{2 \\cdot \\sum_{i < j, \\, i,j \\in C_k} A_{ij}}{|C_k| (|C_k| - 1)}.\n$\nThe inter-cluster density between two groups $C_k$ and $C_l$ ($k \\ne l$) is defined as the proportion of observed edges across the groups:\n$\n\\text{InterDensity}(C_k, C_l) = \\frac{\\sum_{i \\in C_k, \\, j \\in C_l} A_{ij}}{|C_k| \\cdot |C_l|}.\n$This information will be added to the revised supplementary materials.\n\n4. Response to Question 2: \n\nExtensive domain research consistently confirms the inherent sparsity of inter-group connections, a phenomenon substantiated across two primary contexts. In genomic networks (BRCA, LUAD), pathway-to-pathway interactions exhibit natural sparsity since regulatory relationships occur only between specific pathway pairs, such as signaling pathways governing metabolic pathways, while the majority maintain no direct functional interactions [Mariamidze et al. (2018)]. Similarly, within social networks, community-to-community connections demonstrate sparsity due to individuals' predominant engagement within their primary communities, with cross-community interactions remaining comparatively infrequent [Rashidi et al. (2024)]. This foundational assumption of inter-group sparsity aligns with established practices in prior methodological work, such as Cheng et al. (2017) and Shan et al. (2020). Critically, this premise remains consistent with low-level  variable sparsity patterns: when no connections exist between variables across groups, it follows logically that their corresponding group-level interactions should similarly be absent. It should be noted that in scenarios involving minimal group divisions, such as the binary gender classification in Friendship networks, inter-group parameters will not vanish entirely. We will expand upon these considerations in Section 5 of the manuscript.\n\n5. Response to Question 3: \n\nIn our study, we utilized the sample\\_pa function from the R igraph package to stochastically generate networks. This function implements the Barabási-Albert (BA) model, whose core mechanism of preferential attachment represents the classical generative framework for scale-free networks. The resulting networks exhibit node degree distributions that asymptotically follow a power-law distribution. We acknowledge that for small-scale networks, the degree distribution may deviate from a strict power-law due to stochastic fluctuations. Nevertheless, even in small-scale scenarios, the sample\\_pa function remains widely employed for generating randomized modular networks [Ge et al. (2024), Qin et al. (2024)], as it effectively captures two fundamental characteristics of real-world networks: the presence of hub nodes and localized sparsity patterns. To address your comment, the revised manuscript will replace the term \"scale-free\" with the precise descriptor \"Barabási-Albert network\" throughout.\n\nReference:\n\n[1] Cheng, L., et al. (2017). Multilevel gaussian graphical model for multilevel networks. Journal of Statistical Planning and Inference, 190, 1-14.\n\n[2] Codazzi, Laura, et al. (2022). Gaussian graphical modeling for spectrometric data analysis. Computational Statistics \\& Data Analysis, 174, 107416.\n\n[3] Ge, Y., et al. (2024). Structured feature ranking for genomic marker identification accommodating multiple types of networks. Biometrics, 80(4), ujae158.\n\n[4] Lingjærde, C., et al. (2021). Tailored graphical lasso for data integration in gene network reconstruction. BMC Bioinformatics, 22(1), 498.\n\n[5] Liu, H., and Wang, L. (2017). Tiger: A tuning-insensitive approach for optimally estimating gaussian graphical models. Electronic Journal of Statistics, 11(1), 241-294.\n\n[6] Mariamidze, A., et al. (2018). Oncogenic Signaling Pathways in The Cancer Genome Atlas. Cell, 173(2), 321-337.\n\n[7] Navarro, M., et al. (2024). Fair glasso: Estimating fair graphical models with unbiased statistical behavior. Advances in Neural Information Processing Systems, 37, 139589-139620.\n\n[8] Qin, X., et al. (2024). Estimation of multiple networks with common structures in heterogeneous subgroups. Journal of Multivariate Analysis, 202, 105298.\n\n[9] Rashidi, R., et al. (2024). Prediction of influential nodes in social networks based on local communities and users’ reaction information. Scientific Reports, 14(1), 15815.\n\n[10] Shan, L., et al. (2020). Joint estimation of the two-level gaussian graphical models across multiple classes. Journal of Computational and Graphical Statistics, 29(3), 562-579.\n\n[11] Szklarczyk, D., et al. (2025). The STRING database in 2025: protein networks with directionality of regulation. Nucleic Acids Research, 53(D1), D730-D737."}, {"review_id": "Xdwmuh5B9B", "reviewer": "Reviewer_8hqk", "summary": "This paper proposes a deep learning framework for bi-level network inference via node-wise regression for predicting each variable using all others. The authors introduce hierarchical selection layers with adaptive sparsity penalties to capture both group-level and variable-level network structure.", "questions": "See Strengths And Weaknesses section", "limitations": "The limitations discussion is embedded within the conclusion rather than treated separately, and should be expanded.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths: Significance and Originality\n\nThe hierarchical selection mechanism via a group-level layer and a variable-level layer is conceptually appealing and applying neural networks with sparsity penalties to the variable selection setting (particularly to the hierarchical one) appears novel. The theoretical analysis provides formal guarantees and the problem formulation is clear.\n\n## Weaknesses (Quality and Clarity)\n\n- I could not find Stage 1 of the two-stage algorithm from Algorithm 1 in the code. Instead, adaptive weights are computed from random initialization.\n- The experimental evaluation lacks proper validation methodology. Neither the paper nor code describes train/validation/test splits or cross-validation for selecting the super-sensitive sparsity hyperparameters in the baselines, as well as in the proposed method.\n- Baseline choice is hard to interpret: the paper does not discuss its baselines in the related-work section, and the methods mentioned there, i.e., Shan et al. (2020) [22] and Cheng et al. (2017) [4], are not included in the experiments. Comparing against these baselines, or discussing the used baselines in the related work, would strengthen the paper.\n- Placing the related-work section in the appendix is problematic, because clear positioning within the field is central to a scientific paper. The discussion there is quite high-level, mentioning only three competitors without sufficient detail. The authors should include the foundational variable selection work of Meinshausen & Bühlmann (2006), and also cite this when they introduce their variable selection setting, and they should better emphasize their novelty by discussing which aspects have appeared elsewhere. From my (shallow) search, applying neural networks to this variable selection setting appears novel, which would strengthen the paper's contribution if properly positioned in the main text.\n\n### Minor:\nThe architecture uses only 3 layers total, so calling it \"deep\" learning (as in the title) seems overstated. Additionally, referring to a standard linear layer as a \"core innovation\" (line 105) is a bit off.\n\n## Conclusion:\nWhile I appreciate the novel idea and the theoretical contributions look promising, the problems stated above need to be addressed before I can give a positive rating.", "was_revised": true, "final_justification": "Issues Resolved:\n- Code availability: Authors acknowledged the missing Stage 1 implementation and committed to providing it (but missing code is suboptimal of course). Authors will open-source their code upon acceptance.\n- Related work issue: Resolved. Additionally, the code for the closest competitors is not publicly available, which strengthens the authors’ contribution because they will release their own code upon acceptance.\n- Cross-Validation was performed during rebuttal and shows stable results for hyperparameter selection.", "rebuttal": "Dear Reviewer,\n\nWe sincerely appreciate your valuable feedback. We have systematically addressed the three core concerns, including parameter selection mechanisms, method comparison completeness, and paper structure optimization. Below is our point-by-point response:\n\n1. Response to Comment 1:\n\nWe sincerely apologize for the oversight in our code submission, where the implementation corresponding to Stage 1 (initial estimation step) was inadvertently omitted. The current version initializes adaptive weights from random values as a placeholder. In the updated version of our code, which will be made available upon paper revision or acceptance, we will include the full two-stage implementation as described in Algorithm 1. We apologize for the confusion and appreciate your understanding.\n\n2. Response to Comment 2: \n\nGiven the unsupervised nature of our problem, which makes conventional train-validation splitting not applicable for parameter selection via cross-validation, we established the parameter configuration $\\lambda_k=\\zeta_k = c n^{-1/8}$ for $k=1,2$ based on the convergence requirements specified in Theorems 3.6 and 3.7. Here, $c$ denotes a tunable constant and $n$ represents sample size. This formulation aligns with established practices in prior sparsity-inducing models [Liu et al. (2017)]. Computationally, $n^{-1/8}$ yields values between 0.46 and 0.26 for sample sizes ranging from 500 to 50,000. To determine the optimal $c$, we evaluated F1-score performance across values $\\{0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 1, 1.5\\}$ at different sample sizes. The results showed stable performance across a substantial range of $c$ values, with $c = 0.35$ positioned centrally within this stable region. We therefore use $c = 0.35$ in practical implementations to maintain computational efficiency without compromising on accuracy. Supporting results will be presented in the supplementary documentation of the revised manuscript. \n\nRegarding the hyperparameters in our DNN, we conducted exploratory experiments and observed that a configuration of 1000 training epochs, three hidden layers, and 50 nodes per layer offers an optimal balance between computational accuracy and efficiency across both simulated and real-world datasets. Consequently, we adopted this standardized configuration throughout our study. This practice of fixing DNN hyperparameters aligns with established methodologies in prior literature [Sun et al. (2020)], where such standardization enhances reproducibility while maintaining model performance. For the baseline models, we adhered to the parameter selection methodologies recommended in their original publications. We will explicitly document these parameter selection procedures in the supplementary materials. \n\n\n3. Response to Comment 3: \n\nAccording to your suggestions, besides BLGS, we will add the comparative methods Fair Glasso and DeepGRNCS to the Related Works section. Specifically, we will add the reference for Fair Glasso [19] at line 410: \n\n“Notably, Fair Glasso [19] specifically leverages group information to estimate graphical models with provably unbiased statistical behavior, addressing fairness concerns in network inference.”\n\nAdditionally, DeepGRNCS [15] will be added at line 420: \n\n“These methods harness the expressive power of neural networks to model intricate dependencies between variables, often using architectures such as deep neural network with pre-trained model [15], Variational Autoencoders (VAEs) [21] and Graph Convolutional Networks (GCNs) [18].”\n\nFurthermore, although the pioneering bi-level GMM methods (Shan et al., 2020 [22]; Cheng et al., 2017 [4]) mentioned in Appendix A.1 are groundbreaking, their lack of publicly available code makes reproduction difficult. Based on this, this paper selects the latest derivative and improved methods, BGSL and Fair Glasso, as baselines.\n\n4. Response to Comment 4: \n\nDue to space constraints in the main text, we only provided a brief summary of prior work in the first section of the main body (lines 43-54). Consistent with the practice of most papers published at NeurIPS, we have placed a more extensive discussion of related work in Appendix A.1. Following your suggestion, we will expand the description of related work in Appendix A.1 as detailed below: \n\n“To tackle hierarchical network estimation, several GGM-based methods have been developed. Cheng et al. (2017) [4] introduced a multilevel Gaussian graphical model for nested data structures. Shan et al. (2020) [22] proposed a framework for joint estimation of two-level GGMs across multiple classes. Park et al. (2017) [6] focused on learning block-structured graphical models using variable groupings. Notably, Fair Glasso [19] specifically leverages group information to estimate graphical models with provably unbiased statistical behavior, addressing fairness concerns in network inference.”\n\nAdditionally, we will incorporate the work by Meinshausen \\& Bühlmann (2006) you mentioned at line 45 of the main text: \n\n\"Furthermore, a series of methods reformulate the Gaussian Graphical Model (GGM) estimation problem into a set of sparse linear regression models, often offering greater computational efficiency [Meinshausen \\& Bühlmann, 2006].\" \n\nHere, we will specifically highlight its limitation in being able to only capture linear relationships, whereas the proposed method innovatively leverages DNNs to effectively capture nonlinear dependencies among variables.\n\n5. Response to Comment 5: \n\nDiscussion on the number of network layers will be added to the Discussion section. We experimented with deeper network architectures but observed no significant performance improvement despite substantially increased computational costs. Therefore, consistent with existing DNN approaches for high-dimensional omics data and social data [Chen et al. (2021), Lei et al. (2024)], we retain a three-layer architecture. Our model architecture permits seamless extension to deeper networks when application contexts warrant it.\n\nWe sincerely appreciate your valuable suggestions regarding innovation highlights. Compared to standard linear layers, our core innovation lies in leveraging neural networks for variable selection. We will revise the relevant descriptions accordingly and explicitly state the key contribution at the end of the Introduction:\n\n\"Unlike traditional linear methods, this work pioneers the use of neural networks for bi-level variable selection, effectively addressing the challenge of modeling non-linear interdependencies across groups.\"\n\nFollowing your recommendation, we will add a dedicated Section 5.3 \"Limitations\" to discuss the constraints of our approach, thereby enhancing the overall structural clarity of the paper.\n\nReference:\n\n[1] Chen, Y., Gao, Q., Liang, F., \\& Wang, X. (2021). Nonlinear variable selection via deep neural networks. Journal of Computational and Graphical Statistics, 30(2), 484-492.\n\n[2] Lei, Y., Huang, X. T., Guo, X., Hang Katie Chan, K., \\& Gao, L. (2024). DeepGRNCS: deep learning-based framework for jointly inferring gene regulatory networks across cell subpopulations. Briefings in Bioinformatics, 25(4), bbae334.\n\n[3] Liu, H., \\& Wang, L. (2017). Tiger: A tuning-insensitive approach for optimally estimating gaussian graphical models.Electronic Journal of Statistics, 11(1), 241-294.\n\n[4] Sun, T., Wei, Y., Chen, W., \\& Ding, Y. (2020). Genome-wide association study-based deep learning for survival prediction. Statistics in medicine, 39(30), 4605-4620."}, {"review_id": "jVU0aJex88", "reviewer": "Reviewer_snsr", "summary": "In this work, authors propose DeepBLNet, a novel deep learning framework for bi-level network estimation. This method captures nonlinear dependencies and reconstructs inherent hierarchical structures in the data. By integrating hierarchical selection layers and compositional deep architectures, the method identifies both high-level group and low-level variable dependencies while overcoming the limitations of traditional Gaussian and linear models. The paper provides theoretical guarantees for consistent recovery of hierarchical structures and demonstrates superior empirical performance on synthetic and real-world datasets.", "questions": "Here are my questions about the work. Upon receiving satisfying responses to these questions (either in the form of additional analysis/results or explanation as to why the answer isn't important) as well as some assurance relevant information will be included in the updated manuscript, I will revise my score.\n\n1. **Scalability:** How does the method scale with dataset size, specifically in terms of model size? I think it would be valuable to understand how the cost grows for large-scale cases where this work would be most helpful. Authors provide runtimes but do not provide explicit computational complexity or memory requirements.\n\n2. **Interpretability:** Authors claim improved interpretability due to bi-level structure, and provide some insight into how this works in Appendix A3. However, it remains mostly unclear to me how easily domain experts can extract actionable insights from the learned networks. How could one leverage DeepBLNet's so-called interpretability? Can you provide examples for the datasets considered?\n\n3. **Gaussian noise assumption:** How does the assumption of Gaussian noise affect the usefulness theoretical guarantees? Are there relevant examples of datasets/cases where this assumption doesn't hold? How can we understand the value of the method in those cases?", "limitations": "Authors do include a limitations section. I have remaining concerns about model complexity (see Q1) as well as the significance of the Gaussian noise assumption (Q3).", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths**\n- This is a novel, original deep learning approach to bi-level network estimation.\n- The paper provides clear motivation for the work, and specifies the use cases (high dimensional, non-linear) for which it would be most useful.\n- Theoretical guarantees are convincing for the case considered of Gaussian noise. The introduced notion of bi-level consistency appears to be a useful new metric for the field, especially for high-dimensional cases.\n- Empirical results on synthetic and real-world datasets show this is a valuable tool for non-linear, high-dimensional settings. Authors are transparent with the caveat of associated longer runtimes. \n\n**Weaknesses**\n- The approximation theory is rigorous but dense, which is unfortunate because of how useful it could be for practitioners. A more intuitive explanation of the implications (for example, what network depth and width are needed for given data) would improve accessibility.\n- Similarly, the methodology section is not particularly accessible to someone not familiar with network estimation -- Figure 1 is helpful, but would be even more useful if better leveraged in Section 2. As concepts are introduced and steps are outlined, could we follow along in the Figure? Or introduce smaller figures?\n- The work would benefit from additional transparency around model complexity and how that scales (see Q1), practical interpretability (Q2), and effects of the Gaussian noise assumption (Q3).", "was_revised": true, "final_justification": "My main concerns with the work have been largely addressed.\n1. **Scalability:**Authors provide complexity experiments for both runtime and peak memory, that align with theoretical expectations\n\n2. **Interpretability:** Authors provide an example of a dataset (gene pathways) that benefits fromt eh method's so-called interpretability. Authors promise to work this into the paper.\n\n3. **Gaussian noise assumption:** Authors provide a convincing explanation that this assumption is not damaging, and will support this with empirical evidence.\n\nAs such, I raise my score.", "rebuttal": "Dear Reviewer:\n\nSincerely thank you for your profound insights and valuable suggestions on this article. We have systematically addressed the concerns about theoretical assumptions, graphic presentation, computational efficiency, and interpretability. Below is our point-by-point response:\n\n1. Response to Comment 1:\n\nWhile we acknowledge that the expression of the approximation error bound in Theorem 3.2 appears complex, its derivation is built upon the foundational work of Jiao et al. (2023) and has been widely validated in subsequent studies[Jiao et al. (2025), Wang et al. (2024), Shen et al. (2024)]. Crucially, compared to previous bounds exhibiting exponential dependence on dimension $d$ [Lu et al. (2021)], our demonstrates only polynomial dependence on sparsity $s$ as $s^{\\lfloor\\beta\\rfloor+(\\beta\\vee1)/2}$. Furthermore, it achieves the optimal approximation rate $(NM)^{-2\\beta/s}$, explicitly linked to network depth and width (parameterized by $N$ and $M$). This implies that sparser nonlinear functions require smaller network dimensions for accurate approximation, which is an observation aligning with intuition. The bound is fully explicit, containing no hidden or undefined parameters. In addition, in Theorem 3.6 and Theorem 3.7, we establish selection consistency and derive a parametric estimation error bound for the Group Lasso + Lasso estimator and adaptive bi-level sparse estimator respectively. This result features a straightforward expression consistent with prior literature [Dinh et al. (2020)] : the estimation error decays as sample size $n$ increases. These intuitive explanations will be added in the revised paper. \n\n2. Response to Comment 2:\n\nIn order to improve the readability and logicality of the paper, we will migrate Figure 1 to Section 2 as suggested, ensuring that readers can refer to the flowchart immediately when reading the method description, achieving a tight combination of graphics and text for easy understanding. In addition, we will adjust the layout of internal elements in the image, placing the network architecture flowchart in the upper area of the image and the lower area to display the specific network structure. This hierarchical structure is more in line with readers' cognitive habits. On this basis, we will add clear and concise text annotations at each step of the image, and match the formulas with the schematic diagram one by one to ensure that the diagram can convey the core working mechanism of DeepBLNet most intuitively, helping readers to understand the overall framework and key components of DeepBLNet more clearly.\n\n3. Response to Question 1:\n\nWe fully agree on the core importance of computational efficiency in large-scale application scenarios. In order to comprehensively evaluate and clearly demonstrate the computational characteristics of DeepBLNet, we will report the peak memory usage data of the model on existing simulated datasets and multiple real-world datasets in the revised version. In addition, we will add experiments to test the changes in memory usage and running time with sample size $n$. As shown in Table 4 of the paper, the running time shows a significant upward trend with the increase of feature dimension $p$. Based on the structural characteristics of the model, we expect that the time complexity of the model is relatively insensitive to the increase of sample size $n$, and its growth trend is much lower than that with the increase of $p$.\n\n4. Response to Question 2:\n\n\nIn Appendix A3, we mainly explain how to use our method when there are overlaps among different groups. To address your question, in the revised version, we will explain it in two parts through examples. Specifically, we will take the gene networks involved in BRCA and LUAD as examples.\n\nIn the first scenario, assume there are 3 groups, each corresponding to a pathway, indicating that they have similar functions, and the elements in each group do not overlap. When $\\theta_{ll'}$ ($l\\neq l', l,l' = 1,2,3$) is 0, it means that there is no relationship between the $l$-th and $l'$-th pathways. Consequently, the corresponding $j$th and $k$th genes between the groups have $\\theta_{ll'} \\Gamma_{j}^{[:,k]}=0$, so there will be no connections. And assume that $\\theta_{ll}=1$ within the group of the $l$-th pathway, so the connections within the pathway are only determined by $\\Gamma_{j}^{[:,k]}$.\n\nFor the second scenario where, for example, different genes may have multiple functions, the same gene may belong to different pathways and as such, different groups may have overlapping genes. Specifically, let us consider three pathways, denoted as $G_1$, $G_2$, and $G_3$. Each pathway contains five genes: $G_1 = (1,2,3,4,5)$, $G_2 =(4,5,6,7,8)$, $G_3 = (7,8,9,10,11)$. Therefore, $G_1$ and $G_2$ have overlapped genes $4$ and $5$, and $G_2$ and $G_3$ have overlapped genes $7$ and $8$, but $G_1$ and $G_3$ have nothing in common. We can separate the original groups as $G^\\prime_1=(1,2,3), G^\\prime_2=(4,5), G^\\prime_3=(6), G^\\prime_4=(7,8), G^\\prime_1=(9,10,11)$ without overlapping and $G_1=(G^\\prime_1,G^\\prime_2), G_2=(G^\\prime_2,G^\\prime_3,G^\\prime_4),G_3=(G^\\prime_4,G^\\prime_5)$. Denote $\\theta^\\prime_{ll'}=0$ when the reorganized non-overlapping $l$th and $l'$th groups are independent.\n \n Our decision rules for the conditional dependency among the original three groups are defined as follows. For groups with no overlaps such as $G_1$ and $G_3$, if $\\theta^\\prime_{14} = 0$ and  $\\theta^\\prime_{15} = 0$ and  $\\theta^\\prime_{24} = 0$ and $\\theta^\\prime_{25} = 0$, we say $G_1$ and $G_3$ are conditionally independent because all of the variables in $G_1$ are conditionally independent of those in $G_3$. Otherwise, we say $G_1$ and $G_3$ are conditionally dependent.\n For adjacent groups with overlaps, such as $G_1$ and $G_2$, if $\\theta^\\prime_{13} \\neq 0$ or $ \\theta^\\prime_{14} \\neq 0$, we say $G_1$ and $G_2$ are conditionally dependent because at least one variable in $G_1$ unique is conditionally dependent on at least one variable in $G_2$ unique.  If $\\theta^\\prime_{13} = 0$ and $ \\theta^\\prime_{14} = 0$, but $\\theta^\\prime_{12} \\neq 0$ or $\\theta^\\prime_{23} \\neq 0$ or $\\theta^\\prime_{24} \\neq 0$, we say $G_1$ and $G_2$ are conditionally dependent because the common variables are conditionally dependent on at least one variable in both groups. Otherwise, we say $G_1$ and $G_2$ are conditionally independent. In the revised manuscript, we will supplement the above descriptions with a schematic diagram to enhance reader comprehension of the bi-level architecture.\n\n5. Response to Question 3:\n\nThe core role of the normality assumption regarding the error distribution lies primarily in supporting the proof of the Lipschitz property for the empirical risk function in Lemma A.3. This conclusion remains valid under the assumption that the errors follow a sub-Gaussian distribution. Sub-Gaussian distributions constitute a broader class than the normal distribution, encompassing various heavy-tailed distributions. This relaxation significantly enhances the universality of our theoretical results. In the revised appendix, we will provide the complete proof based on the sub-Gaussian assumption. Furthermore, within the Discussion section of the paper, we will add a detailed discussion concerning the necessity of the normality assumption. \n\nConsistent with most existing studies, the theoretical assumptions (e.g., Gaussian/sub-Gaussian distributions) are rarely verifiable in real-world scenarios. The four real-world datasets we employed—Friendship, Co-authorship, BRCA, and LUAD—are widely adopted in related work [Lingjærde et al. (2021), Navarro et al. (2024)], all of which inherently rely on Gaussian assumptions. While our analysis demonstrates relaxation of the Gaussian requirement, the sub-Gaussian assumption remains empirically unverifiable. Nevertheless, comparative evaluation against authoritative ground-truth networks reveals our method achieves decent network reconstruction precision. This provides partial validation of our approach's effectiveness. We will address this empirical evidence in the Discussion section.\n\nReference: \n\n[1] Chen, Y., Gao, Q., Liang, F., \\& Wang, X. (2021). Nonlinear variable selection via deep neural networks. Journal of Computational and Graphical Statistics, 30(2), 484-492.\n\n[2] Dinh, V. C., \\& Ho, L. S. (2020). Consistent feature selection for analytic deep neural networks. Advances in Neural Information Processing Systems, 33, 2420-2431.\n\n[3] Jiao, Y., Shen, G., Lin, Y., \\& Huang, J. (2023). Deep nonparametric regression on approximate manifolds: Nonasymptotic error bounds with polynomial prefactors. The Annals of Statistics, 51(2), 691-716.\n\n[4] Jiao, Y., Kang, L., Liu, J., Lu, X., \\& Yang, J. Z. (2025). Deep approximate policy iteration. The Annals of Statistics, 53(2), 802-821.\n\n[5] Lingjærde, C., Lien, T. G., Borgan, Ø., Bergholtz, H., \\& Glad, I. K. (2021). Tailored graphical lasso for data integration in gene network reconstruction. BMC Bioinformatics, 22(1), 498.\n\n[6] Lu, J., Shen, Z., Yang, H., \\& Zhang, S. (2021). Deep network approximation for smooth functions. SIAM Journal on Mathematical Analysis, 53(5), 5465-5506.\n\n[7] Navarro, M., Rey, S., Buciulea, A., Marques, A. G., \\& Segarra, S. (2024). Fair glasso: Estimating fair graphical models with unbiased statistical behavior. Advances in Neural Information Processing Systems, 37, 139589-139620.\n\n[8] Shen, G., Jiao, Y., Lin, Y., Horowitz, J. L., \\& Huang, J. (2024). Nonparametric estimation of non-crossing quantile regression process with deep requ neural networks. Journal of Machine Learning Research, 25(88), 1-75.\n\n[9] Wang, T., Huang, J., \\& Ma, S. (2024). Penalized Generative Variable Selection. arxiv preprint arxiv:2402.16661."}]

Author rebuttal:

Meta-review: This paper provides a novel extension of Gaussian Graphical models through Deep Networks. It proposes to learn a hierarchical structure through two selection layers(thus the name bilevel) to capture non-linear statistical dependence, beyond the scope of Gaussian Models. Theoretical and experimental results were presented to support the development of the core idea. There was consensus that the paper could be accepted as it adds value to the area of modeling data dependencies through Neural networks.  One of the big concern was code was not made available. If the paper gets accepted the authors should not only incorporate changes discussed during rebuttal and also make code available.
Decision: Accept (poster)