A Latent Multilayer Graphical Model For Complex, Interdependent Systems
Authors: Martin Ondrus
Ivor Cribben
Yang Feng
Venue: NeurIPS 2025 poster ()
Keywords: multilayer networks
neuroimaging
graphical models
latent models
graph learning
Primary area: probabilistic_methods
TL;DR: multiSLICE estimates multilayer networks from multimodal data via a shared latent space, with rigorous theoretical guarantees, extensive simulations, a neuroimaging application, and benchmarks against state-of-the-art methods.

Abstract: Networks have been extensively used and have provided novel insights across a wide variety of research areas. However, many real-world systems are, in fact, a ``network of networks'', or a multilayer network, which interact as components of a larger multimodal system. A major difficulty in this multilayer framework is the estimation of interlayer edges or connections.  In this work, we propose a new estimation method, called multilayer sparse + low-rank inverse covariance estimation (multiSLICE), which estimates the interlayer edges. multiSLICE bridges latent variable Gaussian graphical methods with multilayer networks, offering a flexible framework for modeling processes with irregular sampling and heterogeneous graph structures.  We develop an effective algorithm to compute the estimator. We also establish theoretical conditions for the recoverability of the joint space, analyze how inter-layer interactions influence joint parameter estimation, and provide theoretical bounds on their relationships. Finally, we rigorously evaluate our method on both simulated and multimodal neuroimaging data, demonstrating improvements over state-of-the-art approaches.  Finally, all the relevant R code implementing the method in the article is available on [GitHub](https://github.com/mondrus96/multiSLICE).

A Latent Multilayer Graphical Model for Complex, Interdependent Systems
Martin Ondrus
Ivor Cribben
Neuroscience and Mental Health Institute University of Alberta
Alberta School of Business University of Alberta
Yang Feng
School of Global Public Health
New York University
Abstract
Networks have been extensively used and have provided novel insights across a wide variety of research areas. However, many real-world systems are, in fact, a "network of networks", or multilayer network, which interact as components of a larger multimodal system. A major difficulty in this multilayer framework is the estimation of interlayer edges or connections. this work, we propose new esti- mation method, called multilayer sparse + low-rank inverse covariance estimation (multiSLICE), which estimates the interlayer edges. multiSLICE bridges latent variable Gaussian graphical methods with multilayer networks, offering flexible framework for modeling processes with irregular sampling and heterogeneous graph structures. We develop effective algorithm to compute the estimator. We also establish theoretical conditions for the recoverability of the joint space, analyze how interlayer interactions influence joint parameter estimation, and provide theo- retical bounds on their relationships. Finally, we rigorously evaluate our method on both simulated and multimodal neuroimaging data, demonstrating improvements over state-of-ther art approaches. All the relevant code implementing the method in the article is available on GitHub
Introduction
Our world is complex assemblage of many interdependent systems that can be described as sets of interacting components. Networks (or graphs, which we use interchangcably) succinetly represent these systems, wherein individual actors (nodes) are connected through relationships (edges) that often encode some quantitative meaning (Wasserman and Faust 1994). By representing systems as describe phenomena in social analysis (Wasserman and Faust 1994), genomics (Seal et al. 2023 Argelaguet et al.] 2018). and neuroscience (Bassett and 2017 Bassett et al. 2011). among others.
graphs, local and global graph attributes
Sporns
Gaussian graphical models (GGMs) estimate dependencies such systems by modeling them as multivariate Gaussian distributions, parameterized by the inverse covariance (precision) matrix (Hastie et al. 2209).. However, many real-world systems are in fact a "network networks" (Craven and Wellman 1973) that interact as parts of larger multimodal system. For example, logistics network may contain multiple transportation types or "layers" (c.g., air, rail, pipeline, road), which have within (intralayer) and between (interlayer) relationships. In neuroscience, genes, cells, tissue, anatomical structure, and functional measures can interact and contribute to behaviors or pathologies. Multilayer graphs naturally extend to such cases and encode rich multimodal data into unified representation (Kivelä et al., 2014), but learning the dependence structure remains difficult.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

In such systems, estimating the dependence across different modalities presents challenges due to dif- ferences in the number of measured variables, as well as the sampling rate. neuroscience, for exam- ple, measurement can be made across multiple modalities, such as magntooeletteenepallggrayyy (M/EEG), functional magnetic resonance imaging (fMRI), and structural magnetic resonance imag- ing (sMRI), among others. Each modality has unique strengths and weaknesses that are typically complementary to each other. For example, M/EEG has high temporal resolution but poor spatial resolution, and vice versa for IMRL. Different layers are often tightly coupled, with structure affecting function and vice versa. The challenge combining the information across these modalities that each may differ in the number of variables (nodes) and/or the number of samples. Furthermore, these measurements are taken independently from one another, creating additional issues on how to best structure the "jointness" of these measurements. Ideally, endeavor preserve the natural structure of the data by avoiding up- or down-sampling between different modalities. We may wish to incorporate information from other variables not observed across all modalities. As such, there is clear need for a flexible framework that brings together these multiple modalities.
To this end, we introduce multiSLICE, which bridges multilayer networks (Kivelä et al. 2014) and latent variable GGMs (Chandrasekaran et al. 2012) by providing a flexible framework for modeling processes with irregular sampling and heterogeneous graph structures. In our setting, we consider a sparse + low-rank setup, where each modality has a sparse component, and multiple modalities allowed to exist over a joint low-rank, latent space. Crucially, in a neuroscience context, the sparse component captures local circuits between different regions of the brain, while the low-rank component depicts common environmental influences across modalities (Yatsenko et al., 2015). To our knowledge, this is the first statistical model for estimating interlayer edges in a multilayer network with different numbers of variables and sample sizes, and that is scalable in the number of modalities it can integrate. Our main contributions can be summarized as follows:
We link latent variable Gaussian graphical methods and multilayer networks, establishing a flexible framework to model processes with irregular sampling and non-identical nodes.
We derive an effective algorithm to solve for this estimator.
We establish theoretical conditions for the recoverability of the joint space, analyze how interlayer interactions influence joint parameter estimation, and provide theoretical bounds their relationships.
We rigorously test method both simulated and real experimental data, with compar- isons to the state-of-the-arr, supporting the efficacy of our approach in both settings.
2 Related methods
In general, related methods accommodate either differing numbers of variables (p) between layers or differing sample sizes (n)) Table [ provides summary. Mohan et al. (2014) suggests node- based method (CNJGL) which uses a row-column overlap norm based approach for estimationn Lin et al. (2016) propose multilayer GGM via penalized likelihood estimation (MIGGGM). Gan et al. (2019) uses a Bayesian group regularization method and spike-and-slab I Lasso priors in the proposed BJEMGM method. |Price et al. (202T) suggests cluster fusion regularization (CFR) based method to estimate multiple precision matrices. The BANS method models hicrarchical dependencies using Bayesian node-wise selection (Ha et al., 2021). JMMLE extends this method by decomposing the multilayer problem into two-layer subproblems using neighborhood selection and group-pent pouuppenalized regression (Majumdar and Michailidis 2022). Chang al. (2022) considers graph quilting problem, in which observations of the covariance are missing, and proposes a matrix muarix-completionn based method (LRGQ) estimate the missing entries. Albanese al. (2024) suggests a collaborative graphical lasso (coglasso) for the estimation of multi-omics network data.
CNIGL (2014)
MLGGM (2016)
BJEMGM (2019)
CFR (2021)
BANS (2021)
JMMLE (2022)
LRGQ coglasso (2022)
multiSLICE (Ours)
(2024)
Y2
Yes No
No
Yes
Yes
No
No Yes
Yes
No Yes
Yes
p
Yes
No
No
Yes
Yes
Yes
Table 1: A summary methods in columns with year of publication. Entries indicate whether the approach can handle different sample sizes (n) or different node sets (p) across modalities.

---

3 Preliminaries
3.1 Notation
A graph G is defined as collection of vertices (or nodes) and edges, G (V, For the edge set of weighted graph, we denote the edge between nodes i and j with weight Wiy by the tuple (i.j. way) E E. Altrratively, G can also be described by an adjacency matrix denoted by A. The cardinality a set s, or the number of elements in the set, denoted by |s|. We denote a matrix B by bold uppercase letter, vector b with a bold lowercase letter, and a scalar with a lowercase letter. Bú the ith row and jth column of a matrix B. The vector arising from the jth column B is denoted by bj. Positive definiteness of matrix is denoted by 0 and positive semidefiniteness by > 0. The rank of a matrix B is denoted by RBB. We define a matrix B E RPx as low-rank if R(B) < p. We denote the upper triangle of a square symmetric matrix by (BB), with the corresponding vector by- Pn is a projection operator that selects the indices in Ω2 such that
(Pa(B)lg { Bij if (1,5) EQ if (i,j) Ω
We denote the identity matrix by I, where entries on the main diagonal are and all other entries are O. We denote the singular value decomposition (SVD) of the matrix B by B UAVV AV where the ith value on the diagonal of A is denoted by A; We denote a truncated SVD of rank r by B, = (B) U,A, The norm of a matrix is denoted by -|.
3.2 Multilayer Graph Structure
We define multilayer graph using the principles of Kivelä al. (014)). An l-layer multilayer graph is defined as GM = (Va, EM, V. L), where V is the set of vertices or nodes, L is the set of layers, Vas S x is the set of node-layer tuples. and Es S VMr x VM defines the edge set describing connections between node-layer tuples. A node, u, of a specific layer, a, is encoded as (u, a). The set of intralayer edges, is deñned as FA = {(uu, ar), (u. B)) E FM = 33 while interlayer edges are defined as Ec = {(u, o)) B)) EM a + (} where Ec = EM EA-
supra-adjacency matrix A is a two-dimensional representation of G obtained through a mapping process called "flattening" or "matricization" (Kivelä et al. 2014). For multiSLICE, detailed in Section 3.3] we consider layer-disjoint multilayer network, where each node occupies only one layer. To reduce notation clutter, when a B, we simplify the subscript to a. Hence, an intralayer adjacency matrix is defined as Aa. whereas an interlayer adjacency matrix defined as Aap a + f. Aa,Aa defines the submatrix formed by the intersection of all rows of a and all columns not in a. In general, we use bold subscripts to denote layer-specific quantity.
Remark 1. For a multilayer graph GM with a symmetric supra-adjacency matrix A = AT, we make the following connections to the projection operator Paa
For layer a, the intralayer edges Ea are encoded by projection of that layer, Psta (A), since we have that Ea = {((u, ax), a.) E EEM u, wa).
For I layers, Ω is the union of individual layers' observed indices, Ω = U U (Da) D1 Ul2 U,...,U Q4 The intralayer edges EA across all layers are encoded by Pa(A).
The interlayer edges, Ecc = {(uu, x), ()) E EM u E Wa and u E W3. a + B} are obtained through the complement of Pa. denoted by Poe. Pa=(A) corresponds to the hidden or unobserved edges of A.
The problem of multilayer graph estimation can then be seen as a special case of matrix completion in which A is observed in blocks along the main diagonal. More formally, we observe (a((A) and we wish to recover the full adjacency matrix A by learning the mapping Pa.
To provide some insight, Figure shows a toy example of a 3-layer system with its associated adjacency matrix. In our multilayer system, we are interested in recovering weighted supra- adjacency matrix, in which intralayer edges are organized along the main diagonal in block-like structure and interlayer edges are in the remaining off-diagonal elements.
3

---

Az A2 Azs
Ax Aya Aa
Figure 1: A graphic of weighted adjacency matrix (left), and the associated multilayer graph (right). Intralayer edges shown in the graph solid black lines and color-coded along the main diagonal in block-like structure (A1, A2, A3 in the supra-ade supra-adjacency matrix). Interlayer edges are encoded with dashed lines in the graph (right) and in the matrix (left) the remaining off off-diagonal blocks (uncolored). Note that this matrix is symmetric, which enforces an undirected structure and ensures that interlayer pairs are transposes of each other (A12 A21 T A13 = Asi _ A23 = As2 " The edge between nodes j V in the graph (right) is expressed {(i, (i,2) F.M
3.3 Problem Setting
Assume that L E TRPX is weighted supra-adjacency matrix that composed of square blocks along the main diagonal, corresponding the observed variables. For matrix L", we observe structures related only to intralayer edges and denote the projection into the observed space by Pa- L and its projection (n" are defined as follows:
L;
L
D
L;
L
L
L
0
L
L'
Pa(L')=
Li
Lia
Li
Li
Note that the low-rank latent component corresponding to layer is given by L", Pa, (L") for E (1,...,1), There are also sparse intralayer matrices, the set of which is denoted by {S*)
Ja=1-
We assume that the inverse covariance matrix Σx- originates from the sum of sparse component Sa and a low-rank (or latent, used interchangeably) component Lo such that Σ-- S L where L = Pr. (L") We also assume X X1... X₁ IR XP= are i.i.d. draws from multivariate normal distribution, Xa N(Ha Ex). For the population covariance matrix Ex > and population mean Hoo the finite sample realization is the sample covariance matrix,
1
nh
Σ,
Toi (a)( Ta)
Thx
where Cai and T. are Pax dimensional vectors of the ith sample and sample means, respectively, from the ath layer. Ignoring the l' term, the log-likelihood function is given by
L(S. + La; E.) log det(Sa. + L.) tr tr(E. (Sa + La)).
The joint multilayer likelihood function, then, is given by the sum of individual layerwise likelihoods
Σ C(S.. Pr. EL);
LM(1S) L; {
a=1
By imposing a joint latent matrix L", the number of nodes (Pa) between layers can vary. By focusing on the inverse covariance (precision) matrix rather than the layerwise data vectors (which may be different lengths depending on the number of samples), we do not impose the constraint that the number of samples (na) be identical across layers.
Remark 2.
In short, we model multilayer network as product of layer-specific Gaussian densities whose precision matrices share joint parameterization in the latent space. We set this up by decomposing the

---

precision matrix of each modality into a sparse component unique to that modality and a shared latent component that captures both the within-modalit and the cross-modality edges of the multilayer network. Our objective is to maximize (T) by estimating both the sparse layerwise components as well as the low-rank joint latent component. In Figure we show both the forward data setaegenraating process and the proposed reverse estimation method. The data-generating process comprises two steps: the first step a projection of the latent space into the observed space via P, and the second step is mixing with sparse components to generate the precision matrix. This two-step decomposition important in the next section, where we propose estimating parameters via the reverse process.
s
Pa Pn
L2
2
(S L)
Σ
Σ-
S
Figure An illustration of the data-genenn lata-generating (forward) and parameter estimation (reverse) processes for three-layer system. L denotes the joint latent space parameterized by weighted supra-adjacency matrix, Po denotes the projection into observed layers, and its inverse is Ps The observed S and correspond to edge weights from the intralayer components.
4 Methodology
In (T). the summation over layers decomposes the objective into independent latent variable GGMs. To this end, we introduce the multilayer sparse low-rank inverse covariance estimator (multiSLIEE), which minimizes the following function,
C(Sa: (E=- - Sa
(Sa + La) 1
penalized negative likelihood
covariance fidelity
layers
s.t. RLL) = where <r<p.
Here, p is a tuning parameter for the sparsity in Š and r the pre-specified rank the latent matrix Ż. The objective function is decomposed into a penalized negative log-likelihood term and covariance fidelity term. We propose a two-stage algorithm to solve (2). which follows Figure 2 In the first stage, individual latent variable GGMs estimated, and the second stage, matrix completion is applied to the latent estimates using the block singular value decomposition algorithm (Bishop) and Yu, 2014). Algorithm [anaaamm@eeecce describes an alternating descent algorithm which we refer as sparse low-rank inverse covariance estimation (SLICE). Notice that it can also be used with various other penalties on $ (e.g.. SCAD: Fan et al. 2009) or estimators (e.g.. CLIME; Cai et al. 2011) 1) by simply substituting in a different estimator for S (see the Supplementary Materials for more details).
Algorithm Sparse + low-rank inverse covari- ance estimation (SLICE) with GLASSO
Algorithm Multilayer SLICE (multiSLICE) estimation
Inputs: Σ, p. r, mariter, tol E S 0
Inputs: { p. T, maziter, tol Apply SLICE on E,, for all a
Initialize H =
for i = 1 to maxiter do
i (8) SVD, 16--11
for k to do U, V == SVD, (La)
$ ( GLASSO
(1) ) p)
H. U,A!/2
end for AL AND AS < lol OR ALIS ì) tol Outputs: S, L
end for
= H x H
Outputs: ($). L
5

---

result from Liu We apply SLICE to multilayer networks through coork-coord approach, using it independently for each layer. The estimated L. for each layer is then used for the matrix completion step in the second part of the multiSLICE algorithm. The complete multiSLICE method is detailed in Algorithm 2 For the matrix completion portion, Algorithm first computes per-layer SVDs, L, = A; i == I, , then intersects each pair of latent subspaces by forming Lij U. U. A UJ_ and finally assembles all blocks {La in the global low-rank matrix, L. As an illustrative example, we can consider a two-layer case, where we reconstruct the interlayer latent adjacency matrix L 1 A. U before assembling them into L.
Given assumptions for exact recovery in Theorem 5.1 the rank of each submatrix La must equal that of the overall L' across all a. From Remark we know that depending on Pax and na in each layer, different values of p can be specified for each independent SLICE model. To select and r, we suggest k-fold cross-validation over a grid, where the combined values are based on log- likeliki Simulation experiments elucidating the sensitivity of multiSLICE to these choices are provided in the Supplementary Materials.
Theoretical analysis
We establish consistency conditions for the sparse components and the low-rank joint latent com- ponent for our multiSI.CCE approach. Our proof considers each independent SLICE and matrix completion step separately. Our objective is to show the conditions under which we can exactly recover L" as the sample size each layer, aa -> oo. To recover the full L", we must exactly recover L" from only the observed entries denoted by Pn (L"). ). To do so, ict al. (2017), that describe the general conditions under which exact matrix completion possible. More formally, Liu al. (2017) describe the conditions where Po! exists, which leads us our first result in Theorem 5.1
Theorem 5.1 (Recovery of L"). Let L > with rows and columns indexed by a = {1, p}. Let U U==1 Pa. where Da = wax x Wa and Wa S w. R R(L") Va 1, 1,..., and |u[ [1a=1 Wal [w/, then we have that Pa(L') is invertible.
This result indicates that for exact recovery we only require that each submatrix formed by Pso (L*) be of the same rank as L Next, we require standard assumptions for Sab--aussiantty appropriate regularization, and bounds on the eigenvalues of S" and L' to recover S" and L' (Theorem 5.2).
Remark 3. By Lemma 6.8 Liu et al. (017), the 91/27 -isomeric condition is equivalent to requiring the operators Py- Pa Pu- and Pv- Pa Pv. be invertible. Hence, these invertibility conditions exactly the identifability requirements for the parameter matrix L
Theorem 5.2 (multiSLICE) joint ($} and i consistency). Let L" >> 00 with rows and columns indexed by w = (1,...,p}. Let f = Ua= where Daa = Wa x wax and wa C w. Additionally, let (S 1 1= be the set of true sparse matrices. Then, we have S Saoff Py Siog and L P3 L as na > Ya 1...l
Given the independence of SLICE models, we show the model selection consistency of the sparse component in similar manner. We require standard assumptions for the minimum signal strength of S and irrepresentability. This leads us to the result of Theorem 5.3
Theorem 5.3 (Model Selection Consistency of Sa), For the multiSLICE estimator with LI regular- ization for S in layer a, we have
P (sign(Sa,)) sign( (Sar)) E 5aa) 21 p +
There are further considerations that arise when p and/or n vary between layers.
Remark 4. We apply SLICE models independently, and the conditions for model selection consistency in each layer for the value of p may be different if the number of observed variables in the layer and/or the sample size n different. In such cases, the bounds for ||S - S F, La La lloc are also different. Layers with lower Pa and higher ne therefore require less regularization (through p) and have tighter bounds S and L.

---

Theorem 5.2
establishes the consistency of multiSLICE. However, to establish rate, we require further assumption. In particular, we assume that [|[E Na,al= O La La I)) where for
some layer a, the rate of the unobserved portions the same as the rate of the observed portions, up to constant. With this, we can prove the result in Theorem 5.4
Theorem 5.4 (Rate for joint low-rank, latent space). Let L t 0 with rows and columns indexed by
PP.. Ler f2 = Ua na. where Qa wa x wa and Wa S w. If Ifll L"I Na,allo loe
O(I La La lo)l then we have that
log logPa
1 L'
10
c O==1
for some constant C.
Remark 5. From Theorem 5.4 we have that the oo-norm of the joint space is bounded by the sum of the layerwise Los norms. For two-layer system, the following is implied: The overall bounds on L improve as we increase the sample size from one domain, a 3. even while holding the other number of samples constant. However, we can only recover L exactly when we have exactly recovered La and L
This two-layer intuition can be extended network of any arbitrary size by considering any pair of layers B E I. These phenomena are made evident in the simulation studies in Section 0
6 Simulated data study
1.0
0.3
Simulation 1
-Sore
simel
0.6-
0.0
100
200
800
400
500
100
200
800
400
500
Layers
Simulation 2
F1 Score
sinei u.)
0.8-
0.8
0.1
Rank
1.0
Simulation 3
0.8
|Score
0.8-
1
0.4
0.2
100
200
900
100
300
p.
p.
JMMLE
MLGGM CFR
BANS
BJEMGM
Model
CNUGL
coglasso
LAGQ
multiSi ICE (Ours)
Figure 3: The effects of changing 1, the number of layers (simulation 1) and R(L"), the rank of the latent parameter matrix (simulation 2), and scaling p and the same rate (simulation 3). Each simulation has 100 iterations.
7

---

We consider two challenging simulation studies, which are inspired by structures in neuroimaging (Yatsenko et al. 2015). For each simulation, we generate S. and L". To generate the sparse component for the oth layer, S we first use an initial value for the main diagonal, denoted by B and decay rate, denoted by Ç. The ijth element is defined as
[Beli-al
if Be 2 4, if Be < vb.
S
We then permute S* over rows and columns to randomize the structure. We set 3 = 1.5, ( 2, 0.01, and Pa 100 for all simulations. To generate L", construct p binary matrix Z by, for each row i, selecting one of the 7 columns uniformly at random and setting Zij = 1, with all other
entries in row i equal to zero. Then, we obtain L + 3 x ZZ" For simulation 1, we fix R ) = 2 and vary I; for simulation 2, we fix / and vary R(L") ); for simulation 3, we fix R(L") = 2 and 1 1= 2, and vary p and n jointly. We evaluate estimates of S using the Fl score for sparsity, which typical in the sparse GGM literature (Wang and Allen 2023). ). We evaluate L using the angle between the first eigenvector L and L denoted by sin (( u which is a natural choice for comparing the estimated and true low-rank parameters (Athreya al. 2018).
Ground truth
MLGOM (2016)
BANS (2021)
LRGQ (2022)
muttiSLIEE (ours)
Figure 4: Recovery of L with 2, R(L") = n 500 for methods with non-zero estimates.
Figure 3 shows the results of our simulation study. Across all simulations, we find that multiSLICE outperforms all other methods in the FI score, indicating im- proved estimation of S In addition, it is the only method that converges to- sin ((11, ui) over increasing n.
ward the true L', as indicated by the decrcasing Also, for multiSLICE, the number of layers does not appear to greatly affect the score and sin 0( (ū, uf). while RL') has stronger influence. In Simulation and 2, for all meth- ods, increasing 1 or R(L") has a negative impact on perfor- mance. From the related methods, find that MLGGM has
1000
800
the second-best performance in simulation I, and coglasso has
700
the second-best performance in simulation 2. There is vast
difference in coglasso's s performance between Simulations and 2, suggesting that it is more sensitive to the number of lay- ers than RL) Figure 4 shows the recovered latent parameter matrices.
500
300
200
Figure 5 shows the results of simulating a two-layer system with R[L") 2, with varying 7i and n2, and applying multiSLICE. Our results suggest that increasing the sample size in one layer improves recovery of the joint latent structure n, n2 increases.
100
200 30U
400
5OU CUU 7UD U00 S0U 1003
recovryy
022 0.2
051
Figure 5: L recovery for different layer sample sizes n1. n2 with con- stant 1 2, R(L") =
Lastly, we conduct simulation studies following Wainwright (2009) to validate our theoretical work. In all 100 iterations, we simulate two-layer system, apply multiSLICE, then compute the relevant error norm, and declare recovery successful if that norm falls below constant, predicted by Lemma 9.1. In Figure 6 we plot the empirical success probability against the raw sample size na. and scaled sample size, defined as
nl
Ra
T''s
TíW = Ci log Po
log Pax
8

---

where Sa and Pa are the number of nonzeros and variables in modality a. respectively, and C1 defined in Lemma 9.I We find that by scaling the sample sizes by the appropriate factors, the probability of recovery curves collapse toward common transition point predicted by our theory. Details of hyperparameter selection in all simulations and additional information are in the Supplementary Materials.
S
L
Probability of recovery
0.75
0.90
26
0.25
0.0 5.0
7.5
10.0
0c0
1e:3
2e-3
3e-3
4e-8
6e-3
Figure 6: The probability of successful recovery of S (left panel) and L* (right panel). In each panel, the left subplot x-axis uses the sample size n, whereas the right x-axis uses the scaled sample sizes for S and L", respectively.
Multimodal neuroimaging data study
We apply multiSLICE. and competitor methods multimodal neuroimaging dataset from Wakeman and Henson (2013). In this dataset, 16 subjects are scanned during the presentation of three different facial stimuli. "Famous" and "Unfamiliar" faces are those of people who are publicly well known and those of people who are not, respectively. The "Scrambled" group is a set of images that have the outline or general shape of a face but filled in with white noise; these images serve control stimulus. Each subject has data related to structure (sMRI) and function (MEG, fMRI). Pre- tre-processing steps and further details are the Supplementary Materials, and all experiments are run on a MI MacBook Pro with 16GB of RAM using R 4.4.3.
Each modality has different Pa and fla which is typical in neuroimaging. We exclude methods that assume an equal sample size nix between layers (MLGGM, BANS, JMMLE, coglasso), as this assumption is violated in dataset. Due excessive run-time, CFR omitted. For methods requiring the same Pa across modalities (CNJGL, BJEMGM), we apply SVD to project the data into the smallest joint subspace, corresponding MEG (DMEG = 52). We then apply CNJGL, BJEMGM, LRGQ, multiSLICE to estimate networks independently for each subject and stimulus.
Famous
Unfamiliar
Scrambled
Q($)
H(L)
Q(S)
H(L)
(($)
H(L)
CNJGL (2014)
0,107 (9.65e-03)
N/A
0.106 (1.01c-02)
N/A
0.104 (1.08c-02)
N/A
BJEMGM (2019)
0.084 (1.25c-02)
NIA
0.084 (9.09c-03)
NIA
0.079 (1 (1.27c-02)
N/A
LRGQ (2022)
(1.21e-01)
1.29 (6.63e-02)
0.122 (1.18c-01)
(6.38e-02)
(1.07-001)
1.28 (7.09c-02)
muliSLICE (ours)
0.170 (1.86-03)
0.626 (4.53e-11)
0.171 (1.64e-03) 0.631 (4.02e-02)
0.170 (1.28e-03)
0.660 (1.17e-01)
Table 2: The estimated modularity of the sparse intralayer graphs, Q(S), and the multilayer von
Neumann entropy for the latent supra adjacency matrix, H (L). Bold values are the best results across
methods, highest for ((S) and lowest for H [). for each visual stimulus (Famous, Unfamiliar, and Scrambled). Entries show the mean (standard deviation) across the 16 subjects.
To compare estimates of Š from each method, we use modularity Bullmmore and Sporns 2009). which measures how well communities are separated within the intralayer graphs. This is calculated using Q($) = I (r tr where B = S kk where C a matrix of community assignments, vn is the total number of edges and k is a vector of node degrees. To compare estimates of i from each method, we use the multilayer von Neumann entropy, H, (De Domenico et al. 2015). where lower values are preferred, indicating higher order in the network. The von Neumann entropy is
2m
2m

---

tr (Z) Table 2 shows the results. Across all measures and face stimuli, multiSLICE has the best performance, supporting the efficacy of our method. LRGQ has the second-best performance across all measures, followed by CNJGL and BJEMGM, suggesting that the multiSLICE model more realistic and shows the multiSLICE. adjacency estimates for the different face stimulus conditions. Although they all share low-rank structure, capturing the common task-related variance, each stimulus produces a unique pattern that highlights stimulus- specific multilayer networks.
defined as H(E)= As; In(As,), where Φ is L, normalized to trace of 1, L
yields more favorable estimates of both S and L. Figure shows
Modality
MRI
MEC
sMRI
Famous
Unfamiliar
Scrambled
Figure 7: The estimated supra adjacency matrices, L, obtained from multiSLICE for each visual stimuli presented subject of the Wakeman and Henson (2013) dataset.
8 Discussion
We highlight three key regimes in which our model's guarantees may not hold in practice. The first is that if measurements from different modalities are gathered simultaneously, our assumption of independent modality measurements is violated. Second, any form of temporal or spatial correlation across observations breaches the i.i.d. requirement underpinning our per-modality covariance esti- mators. Lastly, while our sparse, laitt parameterization performs well for neuroimaging data, it remains untested in other application areas. As such, practitioners should validate its suitability before applying to novel domains.
Future studies could adapt multiSLICE to capture sample-to-sampe dependencies via matrix-varis extensions that regularize the row precision matrix in a joint manner. Time- varying and longitudinal domains are also noteworthy, as it is well known that brain imaging data can be non-stationary (Fox al. 2005 Eichele al. 2008: Doucet et al. 2012). For example, for the time-varying case, one could formulate the objective function to jointly estimate {L "T [L]]
with penalties enforcing smooth- ness sparsity in time in Hallac et al. (017) For strictly repeated (longitudinal) measurements, one only needs to treat each time point as a separate "layer" in multiSLICE, ordering covariance inputs by time so that the interlayer estimates {Lr,++1) } capture the within-subjeet evolution directly.
We also consider extending multiLICEE into end-to-end deep-learnng workflows by using its outputs. namely the sparse S and low-rank L estimates define graph inputs for GNNs. For example, Do al. 2023) showed that decoupling graph construction (via Graphical LASSO) from GNN training both accelerates convergence and improves performance, while Sriramulu et al, (2023) 's adaptive dependency learning graph neural networks (ADLNN) uses Graphical LASSO as structural prior before refining edges within GNN. In neuroimaging, Wang et al. (2021); I); Yu et al. (022); y Thapaliya ct al. 2025) utilize graphs representations of regional brain activity for downstream prediction tasks. These sinll-moddlltt methods underscore how graphical models yield effective, interpretable graphs from multivariate data; multiSLICE generalizes this to a multimodal setting, filling a clear gap.
The iterative steps of Algorithm 2 can be unrolled into module and embedded directly into differen- tiable architectures. Alternatively, one can amortize mapping from sample covariances to (S, L) via a small neural network, similar to BBelilossyy et al. (017)), and feed its output into downstream GNNs. In the case of non-Gaussian data, such as raw images, audio, or video, multiSI.CC could be used downstream of modality feature extractors (c.g. CNN).
10

---

Acknowledgements
We would like to sincerely thank the reviewers for their helpful feedback and comments throughout the submission, which have led to several improvements of the paper. The first author was supported by the H. Jean McDiarmaid scholarship provided by the Faculty of Medicine and Dentistry, University of Alberta. The second author was supported by the Natural Sciences and Engineering Research The third author was partially
Council of Canada (NSERC) Discovery Grant RGPIN-2024-06100 supported by the National Science Foundation (NSF) Grant DMS DMS-2324489.
References
Alessio Albanese, Wouter Kohlen, and Pariya Behrouzi. Collaborative graphical lasso. arXiv preprint arXiv:2403.18602, 2024.
Ricard Argelaguet. Britta Velten, Damien Amool, Sascha Dietrich, Thorsten Zenz, John Marioni, Florian Buettner, Wolfgang Huber, and Oliver Stegle. Multi-omics factor analysis framework for unsupervised integration of multi-omics data sets, Molecular Systems Biology, 14(6):c8124) 2018.
Avanti Athreya, Donniell H Fishkind, Minh Tang, Carey E. Priehe, Youngser Park, Joshua T Vogelstein, Keith Levin, Vince Lyzinski, Yichen Qin, and Daniel Sussman. Statistical inference on random dot product graphs: survey. Journal of Machine Learning Research, 18(226): 1-22, 2018.
Danielle S Bassett and Olaf Sporns. Network neuroscience. Nature Neuroscience, 20(3):353-364, 2017.
Danielle S Bassett, Nicholas F Wymbs, Mason A Porter, Peter Mucha, Jean M Carlson, and Scott Grafton. Dynamic reconfiguration of human brain networks during learning. Proceedings of the National Academy of Sciences, 108(18):7641-76 7666766666,, 2011.
Eugene Belilovsky, Kyle Kastner, Gaël Varoquaux, and Matthew B Blaschko. Learning to discover sparse graphical models. In International Conference Machine Learning, pages 440-449 PMLR, 2017.
William E Bishop and Byron M Yu. Deterministic symmetric positive semidefinite matrix completion. Advances in Neural Information Processing Systems, 27, 2014.
Ed Bullmore and Olaf Sporns. Complex brain networks: graph thcoretical analysis of structural and functional systems. Nature Reviews Neuroscience, 10(3):186-198, 2009.
Tony Cai, Weidong Liu, and Xi Luo. A constrained 11 minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106(494):594-607, 2011.
Emmanuel Candes and Benjamin Recht. Exact matrix completion via convex optimization. Commu- nications of the ACM, 55(6): 11-119, 2012.
Venkat Chandrasckaran, Pablo A. Parrilo, and Alan S. Willsky. Latent variable graphical model selection via convex optimization. The Annals of Statistics, 40(4):1935 1967, 2012.
Andersen Chang, Lili Zheng, and Genevera I Allen. Low-rank covariance completion for graph quilting with applications functional connectivity. arXiv preprint arX 2022.
Paul Craven and Barry Wellman. The network city. Sociological Inquiry. 43(3-4):57788, 1973.
Manlio De Domenico, Vincenzo Nicosia, Alexandre Arenas, and Vito Latora. Structural reducibility of multilayer networks. Nature Communications, 6(1):6864, 2015.
Ngoc-Dung Do, Truong Son Hy, and Duy Khuong Nguyen. Sparsity exploitation via discovering graphical models multi-variate time-series forecasting. arXiv preprint arXiv:2306. /7090, 2023.
Gaëlle Doucet, Mikaël Naveau, Laurent Petit, Laure Zago, Fabrice Crivello, Gaël Jobard, Nicolas Delcroix, Emmanuel Mellet, Nathalie Tzourio-Mazoyer Bernard Mazoyer, and Others. Patterns of hemodynamic low-frequency oscillations the brain are modulated by the nature of free thought during rest. Neurolmage, 3994):3194-300, 2012.

---

Tom Eichele, Stefan Debener, Vince D Calhoun. Karsten Specht, Andreas K Engel, Kenneth Hugdahl, D Yves von Cramon, and Markus Ullsperger. Prediction of human errors by maladaptive changes in event-related brain networks. Proceedings of the National Academy of Sciences, 105(16): 6173-6178, 2008
Jianqing Fan, Yang Feng, and Yichao Wu, Network exploration via the adaptive lasso and scad penalties. The Annals of Applied Statistics, 3(2):521, 2009.
Michael D Fox, Abraham Z. Snyder, Justin L. Vincent, Maurizio Corbetta, David C Van Essen, and Mareus E Raichle. The human brain is intrinsically organized into dynamic, anticorrelated functional networks. Proceedings of the National Academy of Sciences, 102((7))9673-98)) 9678, 2005.
Lingrui Gan, Xinming Yang. Naveen Narisetty, and Feng Liang. Bayesian joint estimation of multiple graphical models. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Bue, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
Min Jin Ha, Francesco Claudio Stingo, and Veerabhadran Baladandayuthapani. Bayesian structure learning in multilayered genomic networks. Journal of the American Statisticul Association, 116 (534):605-618, 2021.
David Hallac, Youngsuk Park. Stephen Boyd, and Jure Leskovec. Network inference via the time- varying graphical lasso. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 205-213 2017.
Trevor Hastie, Robert Tibshirani, and Jerome Friedman. The Elements of Statistical Learning: Data mining, Inference, and Prediction, volume 2. Springer, 2009.
Mikko Kivelä, Alex Arenas, Mare Barthelemy, James P Gleeson, Yamir Moreno, and Mason A Porter. Multilayer networks. Journal of Complex Networks, 2(3):203-22) 2014.
Jiahe Lin, Sumanta Basu, Moulinath Banerjee, and George Michailidis. Penalized maximum likelihood estimation of multi-layered gaussian graphical models. Journal of Machine Learning Research, 17(146):1-51, 2016.
Guangean Liu, Qingshan Liu, and Xiaotong Yuan. new theory for matrix completion. Advances in Neural Information Processing Systems, 30, 2017.
Subhabrata Majumdar and George Michailidis. Joint estimation and inference for data integration problems based multiple multi-layered gaussian graphical models. Journal of Machine Learning Research, 23(1): 2022.
Karthik Mohan, Palma London, Maryam Fazel, Daniela Witten, and Su-In Lee. Node-based learning of multiple gaussian graphical models. Journal of Machine Learning Research, 15(1):445448, 2014.
Bradley S Price, Aaron Molstad, and Ben Sherwood. Estimating multiple precision matrices with cluster fusion regularization. Journal of Computational and Graphical Statistics, 30(4):823-834, 2021.
Pradeep Ravikumar, Martin J. Wainwright, Garvesh Raskutti, and Bin Yu. High- nigh-dimenssioal covariance estimation by minimizing LI -enaalll log-determinan divergence. Electronic Journal of Statistics, 5:935 980, 2011.
Adam J. Rothman, Peter J. Bickel, Elizaveta Levina, and Ji Zhu. Sparse permutation invariant covariance estimation. Electronic Journal of Statistics, 2:494 - 515, 2008
Souvik Scal, Qunhua Li, Elle Butler Basner, Laura M Saba, and Katerina Kechris. Refgl: Rapid condition adaptive fused graphical lasso and application to modeling brain region co-expression networks. PLOS Computational Biology, 19(1)ec010758, 2023.
Jakob Seidlitz, František Váša, Maxwell Shinn, Rafael Romero- Comero-Garcia, Kirstie J Whitaker, Petra Vértes, Konrad Wagstyl,, Paul Kirkpatrick Reardon, Liv Clasen, Siyuan Liu, et al. Morphometric similarity networks detect microscale cortical organization and predict inter-individual cognitive variation. Neuron, 97(1):331-24 2018.
12

---

Abishek Sriramulu, Nicolas Fourrier, and Christoph Bergmeir. Adaptive dependency learning graph neural networks. Information Sciences, 625:700-714, 2023.
Bishal Thapaliya, Esra Akbas, Jiayu Chen, Ram Sapkota, Bhaskar Ray, Pranav Suresh, Vince D Calhoun, and Jingyu Liu. Brain networks and intelligence: A graph neural network based approach to resting state fMRI data. Medical Image Analysis, 101:103433, 2025.
Martin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using I1-constrained quadratic programming (lasso). IEEE transactions information theory, 55(5): 2183-2202, 2009.
Martin J Wainwright, John Lafferty, and Pradeep Ravikumar. High-dimensional graphical model selection using I1-regularized logistic regression. Advances in Neural Information Processing Systems, volume 19. MIT Press, 2006.
Daniel G Wakeman and Richard N Henson. A multi-subject, multi-modal human neuroimaging dataset. Scientific Data, 2(1):1-10, 2015.
Lebo Wang, Kaiming Li, and Xiaoping P Hu. Graph convolutional network for fMRI analysis based connectivity neighborhood. Network Neuroscience, 5(1):83-95, 2021.
Minjie Wang and Genevera I Allen. Thresholded graphical lasso adjusts for latent variables.
Biometrika, 110(3)6681-69, 2023.
Stanley Wasserman and Katherine Faust. Social Network Analysis in the Social and Behavioral Sciences, page 3-27. Structural Analysis in the Social Sciences. Cambridge University Press, 1994.
Pan Xu, Jian Ma, and Quanquan Gu. Speeding up latent variable gaussian graphical model estimation via nonconvex optimization. Advances in Neural Information Processing Systems, 30, 2017.
Dimitri Yatsenko, Krešimir Josié, Alexander S Ecker, Emmanouil Froudarakis, R James Cotton, and Andreas S Tolias. Improved estimation and interpretation of correlations in neural circuits. PLOS Computational Biology, 113))=1104088) 2015.
Yue Yu, Xuan Kan, Hejie Cui, Ran Xu, Yujia Zheng, Xiangchen Song, Yanqiao Zhu, Kun Zhang, Razieh Nabi, Ying Guo, al. Learning task-aware effective brain connectivity for fMRI analysis with graph neural networks. In 2022 IEEE International Conference on Big Data (Big Data), pages 4995-4996. IEEE, 2022.
Peng Zhao and Bin Yu. On model selection consistency of lasso. Journal of Machine Learning Research, 7NNov):55412556; 2006.
Tuo Zhao, Han Liu, Kathryn Roeder, John Lafferty, and Larry Wasserman. The huge package for nigh-dimensional undirected graph estimation r. Journal of Machine Learning Research, 13(1): 1069, 2012.
13

---

1. Claims
Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Question:
Answer: [Yes]
Justification:
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 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 attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: We discuss the limitations of our approach in the conclusion. Furthermore, in our theoretical analysis, formal assumptions are outlined, and therefore directly state the constraints of our modelling approach.
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 on 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 a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex. 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 faimess.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
14

---

Justification: We provide assumptions and complete proofs for theoretical results in the Supplementary Materials.
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 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 proofs provided in appendix or 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 of the paper to the extent that it 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: We provide all the steps taken for generating simulations, as well as prepro- cessing the open source neuroimaging data. We provide all code as an Github repository.
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 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 novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way 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 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) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
15

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
Github
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.cm cc/
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 (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips //iiss.cc/ppll wttp:ccpppppccceeeeeeeeee..c.ltt hotol/////gieee/eeeeeuuiiiissss
Coodeooomisssssseeeeeeee www.siiiccciiisiccoooooociiiiiii 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. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary understand the results?
Question:
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level of 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: [Yes]
Justification: For the real data analysis, we provide standard deviation in brackets next to the average values obtained, for each experiment.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for
example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
16

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
OK to report 1-sigma error bars, but one 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 not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures 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: All experiments are run on consumer-grade hardware, and we include timing experiments for our proposed method in the Supplementary Material.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental rums as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (c.g., preliminary or failed experiments that didn't 't make into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., . if there is special consid- eration due to laws or regulations their jurisdiction..
Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification; Our method provides more accurate and robust estimates brain connectivity compared to other oftt--of-the-ar methods, which can accelerate fundamental neuroscience research, improve biomarkers for neurological and psychiatric disorders, and ultimately inform better diagnostics and therapies. However, there are potential negative impacts. For one, there risk of misinterpretation of estimates. Treating the edges correlational network as causal may prompt unsafe interventions. Another concern is privacy. High- resolution connectomes estimated via multilayer networks can, in principle, carry individual- specific signatures. Sharing pooling data without adequate safeguards risks misuse of
17

---

participants' brain data. There are also risks in using this method in unintended ways, such as outside clinical or research contexts (e.g., surveillance of cognitive states). Lastly, there are considerations regarding fairness. If the method is applied to heterogeneous populations without proper care, estimates can systematically misrepresent under-studied groups (e.g.. age, ethnicity). leading to biased conclusions.
Guidelines:
The answer NA means that there is no socictal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
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 to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an 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 that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used 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 unintentional) misuse of the technology.
If there negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms monitor how system learns from feedback over time, improving the efficiency and accessibility of MI...
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have a 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
Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification:
We include appropriate citations for all code and data sets used for the paper.
Guidelines:
The answer NA means that the paper does not existing assets.
18

---

The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license ge.g.. CC-BY 4.0) should be included for each asset. For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
assets are released, the license, copyright information, and terms use in the package should be provided. For popular datasets, papprssitthheeeocc thoode com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-pan 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 assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
We include all code and experiments as Github
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/model part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset 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. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [Yes]
Justification: Although we did not collect any human data part of our neuroimaging experiments, we did use an open source data set which includes human subjects. Details available in the Supplementary Materials.
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 NeurIPS Code of Ethics, workers involved 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 or equivalent for research with human subjects
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?
Question:
Answer: [NA]
19

---

Justification: Although our real data experiments include human subjects, we were not involved in the data collection process all, and are using the open source and de-identified version of the data set.
Guidelines:
The answer NA means that the paper does not involve crowdsoureing nor research with human subjects.
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 NeurlPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods 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, or originality the research, declaration is not required.
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy (ttpg://eeatppsccoooorrrrrreer/p for what should or should not be described.
20


Reviews:
[{"review_id": "DZGBIsST2u", "reviewer": "Reviewer_ioDz", "summary": "This paper introduces multiSLICE, a latent-variable Gaussian graphical model tailored for multilayer networks. The key idea is to decompose each layer’s precision matrix into a sparse intra-layer component plus a shared low-rank component spanning all layers. This low-rank latent factor effectively captures inter-layer connections (influences common across modalities) while the sparse part captures layer-specific associations (e.g. brain-region connections within each modality). Notably, multiSLICE can handle heterogeneous graph structures, without forcing resampling or alignment. The authors propose an efficient two-stage algorithm: first, estimate each layer’s sparse + low-rank structure separately via an EM-like procedure (Algorithm 1), then latent components via a block-wise SVD matrix completion step (Algorithm 2). Finally, the paper empirically validates multiSLICE on both synthetic data (varying number of layers and latent rank) and real multimodal neuroimaging data (MEG + fMRI + sMRI), showing improved accuracy over several state-of-the-art multilayer GGM baselines. Proofs are provided in the supplementary.", "questions": "How might the multiSLICE framework be adapted for time-varying or longitudinal multilayer data?\n\nIt would be insightful to see how multiSLICE performs in a scenario where those methods are applicable. For example, could you run an additional simulation where all layers have equal n and p, and compare multiSLICE with baselines in that controlled setting? This would show that multiSLICE is not only superior in heterogeneous regimes but also competitive in homogeneous ones.\n\nThe method introduces a rank parameter r for L. How should these be chosen in practice? providing guidance or an automatic strategy (even in supplementary) for selecting these hyperparameters would make the method easier to use.", "limitations": "yes", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n-\tCombines sparse + low-rank modeling in a multilayer graphical model formulation.\n-\tAuthors provide a comprehensive theoretical background.\n-\tEmpirical results (synthetic and real neuroimaging data) show superior performance vs. strong baselines.\n\nWeaknesses: \n-\tNotation and matrix completion step could be explained with more intuition and clearly, especially in Algorithm 2, where the authors use block-wise SVD to combine latent components from different layers. It would be helpful to elaborate how this procedure reconstructs a coherent global low-rank matrix from layer-wise estimates and how it addresses issues arising from partial observability or node mismatch across layers.", "was_revised": true, "final_justification": "The response is informative, and it has addressed most of my concerns.\n\nPlease consider to include these additional discussions in the paper.", "rebuttal": "**Weaknesses**\n\n*Notation and matrix completion step could be explained with more intuition and clearly, especially in Algorithm 2, where the authors use block-wise SVD to combine latent components from different layers. It would be helpful to elaborate how this procedure reconstructs a coherent global low-rank matrix from layer-wise estimates and how it addresses issues arising from partial observability or node mismatch across layers.*\n\nThank you for this comment. We assume disjoint node sets across layers, consistent with standard multimodal neuroimaging collection (e.g., no one‑to‑one mapping between EEG channels and fMRI voxels). This motivates the matrix completion step and matches real‑world data, where measurements in modality A and modality B rarely align exactly.\n\nThe matrix completion step in Algorithm 2 follows the block SVD approach from Bishop and Yu (2014). From preliminary simulation studies, comparing block‑SVD against nuclear‑norm penalization and other popular completion methods, we found that this approach achieved the best recovery of missing entries. This also mirrors the results observed in Chang et al. (2022), who found that covariances were more accurately recovered in simulation studies using this method compared to others. To maintain focus, we defer a more exhaustive comparison of alternative approaches to future work.\n\nWe have expanded the explanation of Algorithm 2 to make the intuition more explicit. ``For the matrix completion portion, Algorithm 2 first computes per‑layer SVDs,\n\n$L_i = U_i\\,\\Lambda_i\\,U_i^T,\\quad i=1,\\dots,l,$\n\nthen ``intersects'' each pair of latent subspaces by forming\n\n$L_{ij}=U_i\\,\\Lambda_i^{1/2}\\,\\Lambda_j^{1/2}\\,U_j^T,$\n\nand finally assembles all blocks $\\{L_{ij}\\}$ in the global low‑rank matrix, $L$. As an illustrative case, we can consider a two-layer case, where we reconstruct $L_{12}=U_1\\Lambda_1^{1/2}\\Lambda_2^{1/2}U_2^T$ before assembling them into $L$. This block‑SVD fusion both preserves the latent factors of each layer exactly and aligns their shared structure. As formalized in Theorem 5.1, if the latent matrix of all observed layers has a rank equal to the global rank, then all inter‑layer blocks and hence the full matrix are recovered exactly under our disjoint‑node setup.''\n\n**Questions**\n\n*How might the multiSLICE framework be adapted for time-varying or longitudinal multilayer data?*\n\nThank you for this insightful question. We have added a new discussion in the Conclusion outlining extensions of multiSLICE to time‑varying or longitudinal settings:\n\n``To adapt multiSLICE to a time-varying or longitudinal multilayer data, we suggest the following approaches.\n\nSliding‑windows: Apply multiSLICE over overlapping time windows, yielding a time-varying latent subspace. Although simple and widely used in neuroimaging, this approach may undo some preprocessing steps taken in neuroimaging data analysis (Lindquist, 2024).\n\nChange-point detection: Track a summary statistic of the joint latent factors, such as the principal angles between successive $L$ matrices, and employ offline or online change‑point methods (Cribben et al., 2012; Chen and Zhang, 2015) to segment the time series. Apply multiSLICE separately within each stationary segment.\n\nTemporal regularization: Reformulate the objective function to jointly estimate ${L^{(t)}}_{t=1}^T$, with penalties enforcing smoothness or sparsity in time as in Hallac et al. (2017).\n\nFor strictly repeated (longitudinal) measurements, one only needs to treat each time point as a separate ``layer'' in multiSLICE, ordering the covariance inputs by time so that the inter‑block estimates $\\{L_{t,t+1}\\}$ capture the within‑subject evolution directly.''\n\n*It would be insightful to see how multiSLICE performs in a scenario where those methods are applicable. For example, could you run an additional simulation where all layers have equal n and p, and compare multiSLICE with baselines in that controlled setting? This would show that multiSLICE is not only superior in heterogeneous regimes but also competitive in homogeneous ones.*\n\nThank you for this suggestion. We have added further simulation studies to the Supplementary Materials (Section 12) that compare multiSLICE to the baselines in the controlled setting you have described. In particular, we utilize a two-layer system where $\\mathcal{R}(L^*) = 2$. We then vary $p$ and $n$ equally to compare multiSLICE with BANS, JMMLE, LRGQ, and MLGGM in a homogeneous setting. For brevity, we exclude CNJGL, BJEMGM, CFR, and coglasso, since neither the F1 score nor the $\\sin\\theta$ improved with increasing $n$. The results are outlined in the table below, where each value is the average of 100 replications. We find that multiSLICE, except for the F1 score with $p, n = 50$, has a superior performance over the state-of-the-art methods by a wide margin in this challenging simulation setting.\n\n*Table 1: F1 and* $\\sin\\theta$ *results for various model and sample sizes.*\n\n| $p, n$ | BANS F1    | BANS $\\sin\\theta$ | JMMLE F1   | JMMLE $\\sin\\theta$ | LRGQ F1    | LRGQ $\\sin\\theta$ | MLGGM F1   | MLGGM $\\sin\\theta$ | multiSLICE F1 | multiSLICE $\\sin\\theta$ |\n|:-----|-----------:|------------------:|-----------:|-------------------:|-----------:|------------------:|-----------:|-------------------:|--------------:|------------------------:|\n| 50   | **0.058**  |           0.999   | 0.000      |           0.998   | 0.000      |           1.00    | 0.033      |           0.998   | 0.049         | **0.994**              |\n| 100  | 0.049      |           0.995   | 0.000      |           0.998   | 0.000      |           1.00    | 0.052      |           0.998   | **0.089**     | **0.929**              |\n| 200  | 0.065      |           0.996   | 0.000      |           0.998   | 0.000      |           1.00    | 0.092      |           0.997   | **0.173**     | **0.501**              |\n| 300  | 0.061      |           0.996   | 0.014      |           0.997   | 0.073      |           1.00    | 0.138      |           0.997   | **0.255**     | **0.225**              |\n\n*The method introduces a rank parameter r for L. How should these be chosen in practice? providing guidance or an automatic strategy (even in supplementary) for selecting these hyperparameters would make the method easier to use.*\n\nWe appreciate the reviewer highlighting this point. A full description of our cross-validation procedure for tuning hyperparameters is in the Supplementary Materials (Section 11.1). We have also adjusted the main article, specifically Sections 6 and 7, to better describe our procedure, which is a 3-fold cross-validation over a grid of $\\rho$ and $r$. Hyperparameters are selected as the combination with the best held-out log-likelihood. This is also implemented in our codebase on the anonymous GitHub repository as the `cv.multislice(Xs, folds, rhos, rs)` function, where `Xs`, `folds`, `rhos`, and `rs` are the input data, number of folds for cross-validation, vector of candidate $\\rho$ values, and vector of candidate $r$ values, respectively. A set of simulation studies looking at the effects of $\\rho$ and $r$ on parameter recovery is given in the Supplementary Materials (Section 12).\n\n**References**\n\nWilliam E Bishop and Byron M Yu. Deterministic symmetric positive semidefinite matrix completion. *Advances in Neural Information Processing Systems*, 27, 2014.\n\nAndersen Chang, Lili Zheng, and Genevera I Allen. Low-rank covariance completion for graph quilting with applications to functional connectivity. *arXiv preprint arXiv:2209.08273*, 2022.\n\nHao Chen and Nancy Zhang. Graph-based change-point detection. *The Annals of Statistics*, 43(1):139 – 176, 2015\n\nIvor Cribben, Ragnheidur Haraldsdottir, Lauren Y Atlas, Tor D Wager, and Martin A Lindquist. Dynamic connectivity regression: determining state related changes in brain connectivity. *Neuroimage*, 61(4):907–920, 2012.\n\nDavid Hallac, Youngsuk Park, Stephen Boyd, and Jure Leskovec. Network inference via the time-varying graphical lasso. In *Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining*, pages 205–213, 2017.\n\nMartin A Lindquist. Sliding windows analysis can undo the effects of preprocessing when applied to fmri data. *bioRxiv*, pages 2023–10, 2024."}, {"review_id": "p4MEUL6ips", "reviewer": "Reviewer_JX1i", "summary": "This paper presents multiSLICE, a novel framework for learning multilayer Gaussian graphical models where each layer can have a different number of variables and samples. It models each layer’s precision matrix as a sum of a sparse matrix (intralayer structure) and a shared low-rank matrix (interlayer structure). The two-stage algorithm estimates sparse and latent components per layer, then combines them using a matrix completion approach. The method is theoretically grounded, handles unaligned data, and outperforms prior methods on both simulated and neuroimaging data.", "questions": "Please describe how hyper-parameters are selected", "limitations": "Yes", "rating": 5, "confidence": 1, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n1. Addresses a significant challenge in heterogeneous multilayer graph learning. The method handles differing numbers of variables and samples across layers—an important and previously under-addressed problem. This makes it widely applicable to real-world multimodal settings, such as neuroscience or genomics.\n2. Novel integration of sparse and low-rank modeling with matrix completion. While individually inspired by existing techniques, the combination of latent GGMs, sparse precision estimation, and low-rank matrix completion is original in this context. The model provides both within-layer network structure and explicit interlayer edge estimation.\n3. Scalable, theoretically grounded algorithm. The two-stage procedure (SLICE + matrix assembly) is computationally tractable and supported by recovery guarantees under reasonable assumptions. The design allows each layer to be processed separately, enabling parallelism and better scalability to many layers.\n4. Strong empirical performance across tasks. MultiSLICE consistently outperforms baselines in recovering both sparse and latent structure in simulations. On real neuroimaging data, it produces networks with higher modularity and lower latent entropy, indicating more meaningful and interpretable structure.\n\nWeaknesses:\n\n\nThe paper does not discuss how to choose key parameters like the latent rank or sparsity penalties in practice. This limits immediate usability and may affect performance if the method is misconfigured.\nIncomplete comparisons with recent baselines.", "was_revised": true, "final_justification": "I appreciate the author’s response, which addresses my concern. I continue to uphold my previous rating, leaning toward accepting this paper.", "rebuttal": "**Weaknesses**\n\n*The paper does not discuss how to choose key parameters like the latent rank or sparsity penalties in practice. This limits immediate usability and may affect performance if the method is misconfigured.*\n\nWe appreciate the reviewer highlighting this point. A full description of our cross-validation procedure for tuning hyperparameters is given in the Supplementary Materials (Section 11.1).  We have also adjusted the main article, specifically Sections 6 and 7, to better describe our procedure, which is a 3-fold cross-validation over a grid of $\\rho$ and $r$. Hyperparameters are selected as the combination with the best held-out log-likelihood. This is also implemented in our codebase on the anonymous GitHub repository as the `cv.multislice(Xs, folds, rhos, rs)` function, where `Xs`, `folds`, `rhos`, and `rs` are the input data, number of folds for cross-validation, vector of candidate $\\rho$ values, and vector of candidate $r$ values, respectively. A set of simulation studies looking at the effects of $\\rho$ and $r$ on parameter recovery is given in the  Supplementary Materials (Section 12).\n\n*Incomplete comparisons with recent baselines.*\n\nThank you for this comment. In Section 2, we review eight recent state‑of‑the‑art multilayer and latent‑factor methods, and in Sections 6–7 compare multiSLICE with those methods. Given some of the design differences between methods, such as sample size ($n$) and the number of variables ($p$), we are unable to compare all of them to multiSLICE. If a method does not assume a latent structure, we were unable to compare our latent estimate from multiSLICE. In the real data analysis, because some methods do not posit a latent factor model or require equal sample sizes across layers, direct application to the dataset falls outside their scope. If there are other baselines you think we are missing, we would be happy to include them.\n\n**Questions**\n\n*Please describe how hyper-parameters are selected*\n\nPlease see our response to the first weakness."}, {"review_id": "DIUgONeA4o", "reviewer": "Reviewer_Wcwi", "summary": "This paper addresses the problem of learning the structure of multi-layer network systems by modeling them as Gaussian graphical models, where network connectivity is encoded in the sparsity pattern of the precision (inverse covariance) matrix. The authors formulate the learning task as a regularized maximum likelihood estimation problem, incorporating sparsity-inducing penalties tailored for multi-layer network structures. They derive sample complexity results that establish the number of samples needed for accurate structure recovery and provide high probability bounds on the estimation error of the precision matrix. Additionally, they propose a numerical optimization algorithm to solve the regularized problem efficiently and validate their approach on two benchmark datasets, demonstrating effective structure recovery.", "questions": "1) What is the formal definition of Gaussian distribution on multi-layer networks? If the definition is same as the standard multivariate Gaussian with sparsity pattern of inverse covariance matrix reflecting network structure. Then, the notion of multi-layer network structure is not particularly an interesting probabilistic graphical model. \n2) Is it possible to cast this multi-layer network based graphical model in terms of matrix-variate Gaussian distribution? If so, how? If not why? \n3) Suggestions for improving simulations were already provided in the weakness section. \n4) What is \\mathcal{R}(\\hat{L})=r$ constraint in Eq (2)? Seems like this quantity is not defined in the paper.", "limitations": "The authors, in passing, address as a limitation the fact that their work does not consider the time-varying nature of networks. Personally, I would not consider this as a limitation per se but rather as a natural future avenue. Second, in the checklist, the authors note (paraphrased) that “theorems have assumptions, and those directly state constraints.” I do not completely agree with this viewpoint. Every mathematical theorem has constraints (isn't that the nature of abstract reasoning?). These are not limitations in my view, and hence, I suggest authors to give more serious thought. On societal impact, the authors state there is no potential negative impact. I find this surprising since their simulation work relates to brain studies, where issues such as interpretability, misuse, and privacy often arise.", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 1, "originality": 2, "strengths_and_weaknesses": "Strengths:\nOverall, the paper is clearly written with technically sound mathematical analyses and derivations. While the results align with expectations from existing graphical model literature, the extension to multi-layer networks is well motivated and fairly presented. The results give insights into the interdependency between variables among layers, sparse edges, and the number of samples, thereby further broadening the statistical learning theory of Gaussian graphical models.\n\nWeaknesses:\nThe single biggest weakness is that the contribution feels like a “textbook exercise” to existing theory. I use quotes here to emphasize that while the effort is non-trivial, the motivation for this study feels weak, and the results are expected. I could easily imagine giving this problem as a mini project in a class on non-asymptotic statistics for networks. In the forest of sparsity-regularized likelihoods grown over the last two decades, the current value addition is more in highlighting the need for structure learning and then planting the right tree rather than planting a tree and trying to just the need with very laborious math. However, the motivation itself comes across as purely academic extension, and the final dataset results feel like a throw-off – following the standard template of works in this area.\n\nOn the technical side, since all norms and the cost function are very well-behaved, and the authors resort to fairly standard proof techniques, it is not clear nor well crafted what the theoretical implications really are. Simulation results do not meaningfully (or qualitatively) justify the derived sample complexity statements. It is also unclear if the considered datasets actually adhered to the hypotheses laid out in the theorems. One suggestion is to include a simple but powerful synthetic experiment, similar to what Martin Wainwright did in his 2009 LASSO paper in IEEE IT. Those experiments were simple but spectacularly teased out theoretical results. Once that type of experimental validation is established, the reader gains credibility in your method before seeing brain or cute-cats-running datasets. This is particularly important here because the authors impose their method on certain datasets without any empirical validation that these datasets demanded the proposed setup.", "was_revised": true, "final_justification": "First, the authors did a very detailed and illuminating experiments to support their theoretical results, which was my main concern. Second, they also very clearly addressed the comments of other reviewers. While my initial opinion on the incremental (or textbook exercise) aspect remains the same, nonetheless I've decided to improve the score considering the theoretical analysis, numerical experiments, and the (potential) insights the method suggested for neuroimaging data.", "rebuttal": "**Weaknesses**\n\n*The contribution feels like a \"textbook exercise\" to existing theory, the motivation for this study feels weak, and the results are expected.*\n\nThank you for this comment. We understand that the landscape of sparsity-regularized methodology is saturated. However, our goal is not to add another incremental estimator, but to solve specific scientific problems that current tools leave open.\n\nNeurophysiology studies (e.g., Yatsenko et al., 2015) show that covariances decompose into (i) a *sparse* component, providing undirected associations between brain regions, and (ii) a *low-rank* component, capturing stimuli across modalities. Experimentally, it is of great interest to understand how stimuli influence this low-rank, latent component. Consequently, we adapt this covariance decomposition framework for joint analysis of multiple neuroimaging modalities, and demonstrate its necessity on the Wakeman and Henson (2015) dataset, where existing tools fail to disentangle these structures. We propose the first framework that jointly learns these structures from heterogeneous multimodal data with possibly unequal sample sizes and variable sizes across different data layers. As such, we believe that our method neatly and practically ties together the key areas of multilayer networks, matrix completion, and latent variable Gaussian graphical models.\n\nFurthermore, we provide a key theoretical innovation, a finite‑sample guarantee for joint latent recovery across modalities. In Theorems 5.1 and 5.4, we show that if the latent matrix of each observed modality has the same rank $r$ as the global latent matrix, then all unobserved entries in the full matrix are recovered exactly. Furthermore, we prove that the overall estimation error decomposes cleanly as the sum of the error per layer of each modality. To our knowledge, this is the **first** finite‑sample result for low‑rank completion in a partially‑observed, multilayer setting, effectively bridging the observed and unobserved components of the latent structure.\n\n*Simulation results do not justify the derived sample complexity statements. It is also unclear if the considered datasets actually adhered to the hypotheses laid out in the theorems.*\n\nThank you for this suggestion. We have added further simulation studies to the Supplementary Materials (Section 12), which are inspired by the experiments from Wainwright (2009). We generate two‑layer Gaussian models with a rank of 2 for the latent component, and vary the dimension ($p = 25, 35, 50$) over a range of sample size ($n = 200, 300, 400, 500, 600$). For each setting, we measure the Frobenius‑norm error of the sparse estimate and the infinity‑norm error of the low‑rank estimate, normalizing both by their values at $n = 200$ to isolate the effect of $p$. We find consistency across different $p$ for both sparse and latent components, validating our non‑asymptotic bounds, even in the challenging regime where many latent entries are unobserved. This experiment demonstrates that our theoretical guarantees hold in practice before application to real neuroimaging data.\n\n*Table 2: Mean scaled norms for sparse and latent estimates as a function of variable dimension \\(p\\) and sample size \\(n\\).*\n\n| $$n$$ | $$p=25$$ $$\\|\\hat S - S^*\\|_{F}$$ | $$p=25$$ $$\\|\\hat L - L^*\\|_{\\infty}$$ | $$p=35$$ $$\\|\\hat S - S^*\\|_{F}$$ | $$p=35$$ $$\\|\\hat L - L^*\\|_{\\infty}$$ | $$p=50$$ $$\\|\\hat S - S^*\\|_{F}$$ | $$p=50$$ $$\\|\\hat L - L^*\\|_{\\infty}$$ |\n|:-----|---------------------------------:|--------------------------------------:|---------------------------------:|--------------------------------------:|---------------------------------:|--------------------------------------:|\n| 200  | 1.000                            | 1.000                                 | 1.000                            | 1.000                                 | 1.000                            | 1.000                                 |\n| 300  | 0.882                            | 0.873                                 | 0.885                            | 0.885                                 | 0.888                            | 0.923                                 |\n| 400  | 0.805                            | 0.763                                 | 0.813                            | 0.856                                 | 0.807                            | 0.821                                 |\n| 500  | 0.756                            | 0.748                                 | 0.755                            | 0.785                                 | 0.760                            | 0.744                                 |\n| 600  | 0.703                            | 0.712                                 | 0.702                            | 0.707                                 | 0.713                            | 0.702                                 |\n\n**Questions**\n\n*What is the formal definition of Gaussian distribution on multi-layer networks?*\n\nThank you for this question. We provide a formal definition and derivation of the likelihood in the Supplementary Materials (Section 9.1) which extends the formalism from Section 3.3. We have modified the main article to emphasize these derivations. \n\n\"In short, we model a multilayer network as a product of layer‑specific Gaussian densities whose precision matrices share a joint parameterization in the latent space. We set this up by decomposing the precision matrix of each modality into a sparse component unique to that modality and a shared latent component that captures both the within‑modality and the cross‑modality edges of the multilayer network. This is a more realistic and useful construct, given the experimental evidence cited previously (Yatsenko et al., 2015), and naturally handles heterogeneous $p$ and $n$ across layers.\"\n\n*Is it possible to cast this model as a matrix-variate Gaussian distribution?*\n\nThank you for this insightful suggestion. Casting our model in a matrix‑normal framework is indeed an interesting and promising future direction, as it would allow us to capture both sample‑to‑sample and variable‑to‑variable dependencies. Assuming that each modality follows the matrix-variate Gaussian, \n\n$X_\\alpha \\sim \\mathcal{MN}$  $(0,U_\\alpha,\\Sigma_\\alpha)$ \n\nwe would then impose our sparse + latent decomposition on the column precision as $\\Sigma_\\alpha^{-1} = S_\\alpha + L_\\alpha$. A full matrix normal extension could also regularize the row precision matrix $U_\\alpha^{-1}$ through a similar split — learning the sample‑to‑sample structure — but we leave this rich extension to future work.\n\n*What is $\\mathcal{R}(\\hat{L})=r$ constraint in Eq (2)? Seems like this quantity is not defined in the paper.*\n\nThe operator $\\mathcal R(\\cdot)$ denotes matrix rank.  Hence, $\\mathcal{R}(\\hat{L})=r$ simply enforces that the learned latent matrix $\\hat L$ has rank $r$. The operator first appears and is defined in Section 3.1 - the \"Notation\" portion of \"Preliminaries\". $r$ is also defined, under Equation 2, which reads \"$\\rho$ is a tuning parameter for the sparsity in $\\hat{S}$ and $r$ is the pre-specified rank of the latent matrix $\\hat{L}$.\" In practice, $r$ serves as a hyperparameter controlling the dimension of the low‑rank subspace shared across layers.\n\n**Limitations**\n\n*In the checklist, the authors should give more serious thought to the limitations.*\n\nWe revised the limitations section to:\n\n\"We highlight three key regimes in which our model’s guarantees may not hold in practice.  The first is that if measurements from different modalities are gathered simultaneously, our assumption of independent modality measurements is violated. Second, any form of temporal or spatial correlation across observations breaches the i.i.d. requirement underpinning our per‑modality covariance estimators. Lastly, while our sparse, joint‑latent parameterization performs robustly on neuroimaging data, it remains untested in other application areas; practitioners should validate its suitability before applying it to novel domains.\"\n\n*On societal impact, the authors state there is no potential negative impact.*\n\nThank you for bringing up this concern. We have adjusted the potential impacts statement:\n\n\"Our method provides more accurate and robust estimates of brain connectivity compared to other state-of-the-art methods, which can accelerate fundamental neuroscience research, improve biomarkers for neurological and psychiatric disorders, and ultimately inform better diagnostics and therapies. However, there are potential negative impacts. For one, there is a risk of misinterpretation of estimates. Treating the edges of a correlational network as causal may prompt unsafe interventions. Another concern is privacy. High‐resolution connectomes estimated via multilayer networks can, in principle, carry individual‐specific signatures. Sharing or pooling data without adequate safeguards risks misuse of participants’ brain data. There are also risks in using this method in unintended ways, such as outside clinical or research contexts (e.g., surveillance of cognitive states). Lastly, there are considerations regarding fairness. If the method is applied to heterogeneous populations without proper care, estimates can systematically misrepresent under‑studied groups (e.g., age, ethnicity), leading to biased conclusions.\"\n\n**References**\n\nMartin J Wainwright. Sharp thresholds for high-dimensional and noisy sparsity recovery using l1-constrained quadratic programming (lasso). *IEEE transactions on information theory*, 55(5):2183–2202, 2009.\n\nDaniel G Wakeman and Richard N Henson. A multi-subject, multi-modal human neuroimaging dataset. *Scientific Data*, 2(1):1–10, 2015.\n\nDimitri Yatsenko, Krešimir Josić, Alexander S Ecker, Emmanouil Froudarakis, R James Cotton, and Andreas S Tolias. Improved estimation and interpretation of correlations in neural circuits. *PLoS Computational Biology*, 11(3):e1004083, 2015."}, {"review_id": "9gvfdIcKPt", "reviewer": "Reviewer_pLq3", "summary": "The paper addresses the problem of recovering dependencies among individual systems in multimodal data by modeling them as a multilayer graph, where each system is represented by a latent-variable Gaussian graphical model (GGM). The primary contribution is a statistical model that estimates interlayer edges across disjoint GGMs that vary in structure, number of nodes, and sample sizes.\n\nThe authors provide theoretical guarantees for the recoverability of the joint space, exploring how interlayer interactions impact joint parameter estimation. They also derive theoretical bounds on these relationships. To the best of their knowledge, this is the first model capable of estimating interlayer connections in a multilayer network that handles heterogeneous dimensions, sample sizes, and scales to multiple modalities.", "questions": "N/A", "limitations": "Yes. The limitations have been discussed.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "### **Strengths:**\n\nThe paper is well written and well motivated, with a clear problem formulation. The proposed framework is flexible and shows superior empirical performance compared to baseline methods with insightful ablation studies. It includes theoretical recovery guarantees for the joint parameters. The approach can also handle heterogeneous systems and is scalable to multiple modalities.\n\n### **Weaknesses:**\n\nThe paper makes several assumptions which may limit generalizability in real-world data. These assumptions include: \n1. the individual models are disjoint Gaussians \n2. known intralayer connectivity and assumes that data is generated from a self-sufficient model with that known structure. From my understanding from Line 138, the data is assumed to be not influenced by the interlayer latent dependencies. \n3. prior knowledge of the maximum number of variables (p), including latent ones\n4. decomposability of the joint likelihood into layer-wise summations for tractability. The decomposition is made with few justifications. It would have been more useful if the authors could provide more convincing justifications for this design choice if there were any.\n\nWhile the paper provides a solid solution to the problem at hand, in all honestly, I am unsure whether it fits into the landscape of NeurIPS since the authors solely employ traditional statistical methods and the paper lacks application of neural networks. The authors should enhance the paper's relevancy and impact by demonstrating applications in deep learning contexts to align with contemporary research directions.", "was_revised": true, "final_justification": "Following my discussion with the authors, I decide to improve the rating to Borderline Accept given their additional discussions on the applicability of the proposed method in deep learning frameworks. The borderline evaluation pertains to the limit of the current discussion to high-level ideas, rather than a detailed technical discussion that should have been included in the submitted version. Overall, the paper itself is a sound statistical work itself. Although the contribution to me remains an incremental contribution to the statistics community, I think the problem setup and the completeness of the results, especially regarding the recovery of latent factors, can somewhat inspire the design of new learning algorithms.", "rebuttal": "**Weaknesses**\n\n*The paper makes several assumptions which may limit generalizability in real-world data. These assumptions include: The individual models are disjoint Gaussians with known intralayer connectivity and assumes that data is generated from a self-sufficient model with that known structure. From my understanding from Line 138, the data is assumed to be not influenced by the interlayer latent dependencies.*\n\nThank you for this comment. We emphasize that in multiSLICE, no layer is treated in isolation — each layer utilizes the same joint low‑rank latent matrix, hence its precision matrix combines that shared block with a layer‑specific sparse component. As a result, the structure and parameters estimated in any one layer are directly influenced by the information from all other layers via the common latent structure.\n\nAlthough we focus on the Gaussian case for clarity, multiSLICE naturally extends to non‑Gaussian settings. For example, we could use the approach of Liu et al. (2009) and extend our method to a semi-parametric domain. Further extensions to exponential family graphical models are also worth exploring (Yang et al., 2015).. We have added a note about these potential extensions to the discussion.\n\n*Prior knowledge of the maximum number of variables (p), including latent ones*\n\nFor multiSLICE, no extra prior knowledge of $p$ or the latent dimension is required beyond the data itself. The total variable count $p$ is simply the sum of measured features across modalities. The latent rank $r$ is treated as a tunable hyperparameter rather than an assumed constant: we select $r$ (and the sparsity penalty $\\rho$) via 3‑fold cross‑validation on held‑out likelihood. \n\nIn response to your comment, we have adjusted the main article, specifically Sections 6 and 7, to better describe the procedure we use.  This is also implemented in our codebase on the anonymous GitHub repository as the `cv.multislice(Xs, folds, rhos, rs)` function, where `Xs`, `folds`, `rhos`, and `rs` are the input data, number of folds for cross-validation, vector of candidate $\\rho$ values, and vector of candidate $r$ values, respectively. A set of simulation studies looking at the effects of $\\rho$ and $r$ on parameter recovery is given in the Supplementary Materials (Section 12).\n\n*Decomposability of the joint likelihood into layer-wise summations for tractability. The decomposition is made with few justifications. It would have been more useful if the authors could provide more convincing justifications for this design choice if there were any.*\n\nThank you for this comment. This design choice is not merely for tractability; it is a core feature motivated by established findings in neurophysiology. As shown by (Yatsenko et al., 2015) and others, neural population covariance naturally decomposes into a sparse, neuron-specific component and a shared, low-rank component modulated by brain state. Our likelihood decomposition is a direct mathematical implication of this principle. Specifically, we provide a formal definition and derivation of likelihood in the Supplementary Materials (Section 9.1), which extends the description from Section 3.3. We have adjusted the main article to emphasize these derivations. \n\n\"In short, we model a multilayer network as a product of layer‑specific Gaussian densities whose precision matrices share a joint parameterization in the latent space. We set this up by decomposing the precision matrix of each modality into a sparse component unique to that modality and a shared latent component that captures both the within‑modality and the cross‑modality edges of the multilayer network. This is a more realistic and useful construct, given the experimental evidence cited previously (Yatsenko et al., 2015), and naturally handles heterogeneous $p$ and $n$ across layers.\"\n\n*While the paper provides a solid solution to the problem at hand, in all honestly, I am unsure whether it fits into the landscape of NeurIPS since the authors solely employ traditional statistical methods and the paper lacks application of neural networks. The authors should enhance the paper's relevancy and impact by demonstrating applications in deep learning contexts to align with contemporary research directions.*\n\nThank you for raising this point. We respectfully argue that foundational statistical methodology, particularly for graphical models, remains a vibrant and essential part of the NeurIPS landscape. This is reflected in a significant number of recent papers at NeurIPS (Zhang et al., 2024; Negri et al., 2023; Navarro et al., 2024) and many in the broader literature (Roverato and Nguyen, 2024; Tsai et al., 2022; Wang and Allen, 2023; Chang et al., 2022; Kundu et al., 2019).\n\nIn addition, our work addresses a critical need in hypothesis-driven scientific fields like neuroimaging, where high-dimensional data and limited sample sizes often make deep learning approaches impractical or prone to overfitting (Wodeyar and Srinivasan, 2022; Smith et al., 2011).  In these data-scarce regimes, theoretically grounded statistical models like multiSLICE are indispensable for generating reliable and interpretable insights.\n\nFar from being disconnected from contemporary research, we see this work as enabling future deep learning applications. By providing a robust method for learning structured, low-dimensional representations from complex multimodal data, our work can provide the basis for more sophisticated deep learning architectures.  We believe this contribution aligns well with NeurIPS's mission of advancing machine-learning methodology, from foundational theory to practical application.\n\n**References**\n\nAndersen Chang, Lili Zheng, and Genevera I Allen. Low-rank covariance completion for graph quilting with applications to functional connectivity. *arXiv preprint arXiv:2209.08273*, 2022.\n\nSuprateek Kundu, Joshua Lukemire, Yikai Wang, and Ying Guo. A novel joint brain network analysis using longitudinal alzheimer’s disease data. *Scientific reports*, 9(1):19589, 2019.\n\nHan Liu, John Lafferty, and Larry Wasserman. The nonparanormal: semiparametric estimation of high dimensional undirected graphs. *Journal of Machine Learning Research*, 10(10), 2009.\n\nMadeline Navarro, Samuel Rey, Andrei Buciulea, Antonio G Marques, and Santiago Segarra. Fair glasso: Estimating fair graphical models with unbiased statistical behavior. *Advances in Neural Information Processing Systems*, 37: 139589–139620, 2024.\n\nMarcello Massimo Negri, Fabricio Arend Torres, and Volker Roth. Conditional matrix flows for gaussian graphical models. Advances in *Neural Information Processing Systems*, 36:25095–25111, 2023.\n\nAlberto Roverato and Dung Ngoc Nguyen. Exploration of the search space of gaussian graphical models for paired data. *Journal of Machine Learning Research*, 25(92):1–41, 2024.\n\nStephen M Smith, Karla L Miller, Gholamreza Salimi-Khorshidi, Matthew Webster, Christian F Beckmann, Thomas E Nichols, Joseph D Ramsey, and Mark W Woolrich. Network modelling methods for fmri. *Neuroimage*, 54(2): 875–891, 2011\n\nKatherine Tsai, Oluwasanmi Koyejo, and Mladen Kolar. Joint gaussian graphical model estimation: A survey. *Wiley Interdisciplinary Reviews: Computational Statistics*, 14(6):e1582, 2022.\n\nMinjie Wang and Genevera I Allen. Thresholded graphical lasso adjusts for latent variables. *Biometrika*, 110(3):681–697, 2023.\n\nAnirudh Wodeyar and Ramesh Srinivasan. Structural connectome constrained graphical lasso for meg partial coherence. *Network Neuroscience*, 6(4):1219– 1242, 2022.\n\nEunho Yang, Pradeep Ravikumar, Genevera I Allen, and Zhandong Liu. Graphical models via univariate exponential family distributions. *The Journal of Machine Learning Research*, 16(1):3813–3847, 2015.\n\nDimitri Yatsenko, Krešimir Josić, Alexander S Ecker, Emmanouil Froudarakis, R James Cotton, and Andreas S Tolias. Improved estimation and interpretation of correlations in neural circuits. *PLoS Computational Biology*, 11(3):e1004083, 2015.\n\nShengjun Zhang, Xin Fei, Fangfu Liu, Haixu Song, and Yueqi Duan. Gaussian graph network: Learning efficient and generalizable gaussian representations from multi-view images. *Advances in Neural Information Processing Systems*, 37:50361–50380, 2024."}]

Author rebuttal:

Meta-review: The paper proposes multiSLICE, a latent-variable Gaussian graphical model for multilayer networks, combining sparse within-layer structure with a shared low-rank latent component.

Reviewers’ main concerns about clarity, theory - experiment alignment, and practical usability have been addressed.

For the camera-ready, the authors should (1) integrate the added simulation results and probability-of-recovery plots, (2) include the identifiability discussion explicitly, (3) expand the guidance on hyperparameter tuning and practical usage, and (4) add the discussion on links to deep learning and broader applications.
Decision: Accept (poster)