Learning Graph Structure from Convolutional Mixtures
Authors: Max Wasserman
Saurabh Sihag
Gonzalo Mateos
Alejandro Ribeiro
Venue: Accepted by TMLR ()

Abstract: Machine learning frameworks such as graph neural networks typically rely on a given, fixed graph to exploit relational inductive biases and thus effectively learn from network data. However, when said graphs are (partially) unobserved, noisy, or dynamic, the problem of inferring graph structure from data becomes relevant. In this paper, we postulate a graph convolutional relationship between the observed and latent graphs, and formulate the graph structure learning task as a network inverse (deconvolution) problem. In lieu of eigendecomposition-based spectral methods or iterative optimization solutions, we unroll and truncate proximal gradient iterations to arrive at a parameterized neural network architecture that we call a Graph Deconvolution Network (GDN). GDNs can learn a distribution of graphs in a supervised fashion, perform link prediction or edge-weight regression tasks by adapting the loss function, and they are inherently inductive as well as node permutation equivariant. We corroborate GDN's superior graph learning performance and its generalization to larger graphs using synthetic data in supervised settings. Moreover, we demonstrate the robustness and representation power of GDNs on real world neuroimaging and social network datasets.

Published in Transactions on Machine Learning Research (06/2023)
Learning Graph Structure from Convolutional Mixtures
Max Wasserman Deportment of Computer Science University of Rochester
mwwsserf@ccs.rccesseer.c
Saurabh Sihag
sitaggs@penmmmddciinn..pppp upenn. edu
Department of Electricol and Systems Engineering University of Pernsyloania
Gonzalo Mateos
smateosh@ecc.rrchhesee
Deportment of Electricol and Computer Engineering University of Rochester
Alejandro Ribciro
aribesiro Aseas.upenn. edu
Department of Electrical and Systems Engineering University of Pennsyleania
Reviewed on OpenReview: https: // opt.reveeerrvv.....mmmmm net/ forum? id= OILbPOWERR
Abstract
Machine learning frameworks such as graph neural networks typically rely on given, fixed graph to exploit relational inductive biases and thus effectively learn from network data. However, when said graphs are (partially) unobserved, noisy, dynamic, the problem of inferring graph structure from data becomes relevant. In this paper, we postulate graph convolutional relationship between the observed and latent graphs, and formulate the graph structure learning task as a network inverse (deconvolution) problem. In lieu of aggnddcomposition-bssee spectral methods or iterative optimization solutions, we unroll and truncate proximal gradient iterations to arrive a parameterized neural network architecture that we call Graph Deconvolution Network (GDN). GDNs learn distribution of graphs in supervised fashion, perform link prediction or edge weight regression tasks by adapting the loss function, and they are inherently inductive as well as node permutation equivariant. We corroborate GDN's superior graph learning performance and its generalization to larger graphs using synthetic data in supervised settings. Moreover, we demonstrate the robustness and representation power of GDNs on real world neuroimaging and social network datasets.
1 Introduction
Inferring graphs from data to uncover latent complex information structure is a timely challenge for geometric
deep learning (Bronstein et al. 2017 and graph signal processing (GSP) (Ortega et al. 2018 But it also an opportunity, since network topology inference advances (Dong et al. 2019 Mateos et al. 2019 Giannakis et al. 2018)) could facilitate adoption of graph neural networks (GNNs) even when no input graph available, e.g. for discase prediction in healtheare settings (aaii et al. 2023 or, 3D point cloud segmentation in computer graphies (Wang al. 2019). The problem is also relevant when given graph is too noisy perturbed beyond what stable (possibly pre-trained) GNN architectures can effectively handle (Gama et al. 2020a). Early empirical evidence suggests that even when a graph is available, the structure could be further optimized for a downstream task e.g., when parameterizing a GNN used for node classification in citation, social, and protein-protein interaction networks (Zhao et al.) 2021). Sometimes the observed graph's 's large odge set can make downstream task computationally expensive, which motivates seeking a sparser graph which retains properties of the original one thus making the task feasible (Spielman & Srivastava 2011

---

Published in Transactions on Machine Learning Research (06/2023)
In this paper, we posit convolutional model relating observed and latent undirected graphs and formulate the graph structure learning task as a supervised network inverse problem; see Section 2 for a formal problem statement. This fairly general model is motivated by various practical domains outlined in Section 3 such as identifying the structure of network diffusion processcs (Segarra et al., 2017: Pasdcloup et al.) 2018). as
well as network deconvolution and denoising (Feizi et al. 2013). We propose a parameterized neural net (NV) model, termed graph deconvolution network (GDN), which we train in a supervised fashion to learn the distribution of latent graphs. The node permutation equivariant architect is derived from the principle of algorithm unrolling used to learn fast approximate solutions to inverse problems (Gregor & LeCun 2010 Sprechmann et al. 2015 Monga et al., 2021). an idea. that is yet to be fully explored for graph structure identification. Since layers directly operate on, combine, and refine graph objects (instead of nodal features), GDNs are inherently inductive and can generalize to graphs of different size. This allows the transfer of learning on small graphs to unseen larger graphs, which has significant implications in domains like social networks and molccular biology (Ychudai ct al., [021]. Our experiments demonstrate that GDNs are versatile to accommodate link prediction or edge-weight regression aspects of learning the graph structure, and achieve superior performance over various competing alternatives. Building on recent models of functional activity in the brain a diffusion over the underlying anatomical pathways (Abdelnour et al. 2014 Liang & Wang 2017), we show the applicability of GDNs to infer brain structural connectivity from functional networks obtained from the Human Connectome Project-Young Adult (HCP-YA) dataset. We also use GDNs to predict Facebook ties from user co-location data, outperforming relevant baselines.
Related work. Graph structure learning (a.k.a. network topology inference) has a long history in statis- tics Dempster 1972). with noteworthy contributions for probabilistic graphical model selection; see e.g. (Ko- laczyk] 2009 Friedman et. al. 2008 Drton & Maathuis 2017). Recent advances were propelled by GSP
insights through the lens of signal representation (Dong et al.
2019
Mateos et al.) 2019 Giannakis al.,
2018], exploiting models of network diffusion (Daneshmand et al. 2014] 4) or else leveraging cardinal properties
of network data such smoothness (Dong et al.) 2016 Kalofolias 2016) and graph stationarity (Segarra)
et al. 2017: Pasdeloup et al.) 2018). These works formulate (convex) optimization problems one has to solve for different graphs, and can lack robustness to signal model misspecifications. Network deconvolution approaches in (Segarra et al.
2017 Feizi et al. 2013) operate in the graph spectral domain, and may face scalability issues because they rely on computationally eexptaieii eigendecompositions of the input graph for each problem instance. Morcover, none of these methods advocate a superviscd learning paradigm using NNs to predict adjacency matrices as we propose here. In the broader context of NN-based models, so-termed latent graph learning has been shown effective in obtaining better task-driven representations of relational data for machine learning (ML) applications (Wang et al. 2019 Kazi et al. 2023 Veličković et al.) 2020), or to learn interactions among coupled dynamical systems (Kipf et al.) 2018). With regards to NN architectural design, algorithm unrolling has only recently been adopted for graph structure learning in (Shrivastava et al.) 2020 Pu et al.) 2021); but for different and more specific problems subsumed by the network deconvolution framework dealt with here. Indoed, (Shrivastava et al., 2020) focuses on Gaussian graphical model selection in supervised learning setting and (Pu et al.) 2021 learns graph topologies under a smoothness signal prior.
Summary of main contributions. In this work, we introduce the following graph MT. innovations:
We rntrod GDNs, supervised learning NN model capable of recovering sparse latent graph structure from observations of its convolutional mixtures, ie., related graphs containing spurious, indirect relationships.
The novel supervised setting for graph structure learning is motivated via several real-world applications.
. The node permutation equivariant GDN architecture is derived from the principle of algorithm unrolling to learn fast approximate solutions to inverse problems. The unrolling offers explicit control on complexity (leading to fast inference times) and can seamlessly integrate domain-speeific prior information about the unknown graph distribution. In constructing GDNs we leverage Multi-Input Multi-Output (MIMO) filters- as used in convolutional (C)NNs for further expressiveness and training stability.
On synthetic data, GDNs outperform comparable methods on link prediction and edge- regression tasks across different random-graph ensembles, while incurring a markedly lower (post-training) computational cost and inference time. GDNs are inductive and learnt models transfer across graphs of different sizes. We empirically verify they exhibit minimal performance degradation even when tested on graphs 30x larger.
2

---

Published in Transactions on Machine Learning Research (06/2023)
Finally, using GDNs we propose novel ML pipeline to learn whole brain structural connectivity (SC) from functional connectivity (FC), challenging and timely problem in network neuroscience that motivates the supervised learning setting advocated here. Results on the HCP-YA imaging dataset show that GDNs perform well on specific brain subnetworks that are known to be relatively less correlated with the corresponding FC due ageing-related effects a lestament to the model's robustness and expressive power. Overall, results here support the promising prospect of using graph ML to integrate brain structure and function.
2 Problem Formulation
In this work we study the following nelwork inverse problem involving undirected and weighted graphs G(V,E), where V= NN the set of nodes (henceforth common to all graphs), and C x V collects the edges. We get to observe a graph with symmetric adjacency matrix Ao E RNXN, that is related to a latent sparse, graph Ag R RN of interest via the forward data model
Ao = L1 - h1 Ar + h2A + A
(1)
indirect relationships
for some S N 1 by the Goyley-Hamilton theorem. Matrix polynomials (A; == Ex = h4 A with coefficients h [ho....: Nx]T RK+1 as in the right-hand-side of eq. (1), are known shift-invariant graph convolutional filters; see e.g.: (Orlega et al. 2018 Gama et al.) 2020bb. We postulate that Ao = H (Ar; h) for some filter order K and its associated coefficients h, such that we can think of the observed network as generated via a graph convolutional process acting on AL. That is, one can think of Ao as a graph containing spurious, indirect connections generated by the terms including higher-order powers of latent graph A - the graph of fundamental relationships. More pragmatically Ao may correspond a noisy observation of H(A;; h), and this will be clear from the context when e.g,, we estimate Ao from data.
Recovery of the latent graph Ay is challenging endeavour since we do not know HAL;;)), namely the parameters K or h; see Appendix A.T for issues of model identifiability and their relevance to the problem dealt with here. Suppose that AL is drawn from some distribution of sparse graphs, e.g., random geometric graphs or structural brain networks from a relatively homogencous dataset. Then given independent training samples T= = (46) 1T IT adhering to eq. 1 our goal is to learn parametric mapping & that predicts the graph adjacency matrix Àg (Aoo +) by minimizing a loss function
I( = T
1 Σ
((A 1 ST
(2)
to search for the best parameters O or mean-squared error (MSE) for the more challenging edge-weight regression problem; see Appendix A.5 Notice that eq. | postulates a common filter across graph pairs in T a simplifying assumption shown be tenable in functional magnetic resonance imaging (fMRI) studies [Abdelnour et al.) (2144) where subject-level coupling between the structural and functional modalities (SC-FC coupling) exhibit limited variability. This
The loss é is adapted to the task at hand hinge loss for link prediction
will be relaxed when customize the GDN architecture to learn multiple filters; see Section 4.3 for details.
Potential application domains for which supervised network data T ( ( 1T is available include bioinformaties (infer protein contact structure A, from mutual information graphs Ao of the covariation of amino acid residues (Fezii et al. 2013)) ), gene regulatory network inference from microarray data (see e.g. the DREAM5 project to reconstruct networks for the E. coli bacterium and single-cell eukaryotes), social and information networks (e.g,, learn to sparsify graphs (Spielman & Srivastava) 2011 to unveil the most relevant collaborations A in a social network encoding co-authorship information Ao [Segarra et al. 2017 )), and cpidemiology (such as contact tracing by deconvolving the graphs that model observed discase spread population). In Section 5.2 we experiment with social networks and the network neuroscience problem described next.
3

---

Published in Transactions on Machine Learning Research (06/2023)
3 Motivating Application Domains and the Supervised Setting
Here we outline several graph structure learning tasks that can be cast as the network inverse problem in eq- I and give numerous examples supporting the relevance of the novel supervised setting.
Graph structure identification from diffused signals. Our initial focus is on identifying graphs that explain the structure of class of network diffusion processes. Formally, let æ E IN be a graph signal (i.e., a vector of nodal features) supported on a latent graph & with adjacency AL. Further, let w be a zero-mean white signal with covariance matrix En E Ewww = I. We say that A represents the structure of the signal x if there exists linear network diffusion process in g that generates the signal æ from w, namely x= = >i= A4 w = (Ar; h)w. This is a fairly common generalive model for random ner processes (Barrat et al., 2008 DeGroot 1974). We think of the edges of G as direct (one-hop) relations between the elements of the signal æ. The diffusion described by (A,; h) generates indirect relations. In this context, the graph structure learning problem is to recover a sparse AL from a set X (2;);= of P samples of æ (Segarra et al. 2017). Interestingly, from the model for x it follows that the signal covariance matrix Σ,- = æ = H (Az; E www" 1| (A);)) = H² h)) is also polynomial in A precisely the relationship prescribed by eq. [1 when Ao = Ez. In practice, given the signals in X' one would estimate the covariance malrix, e-g. via the sample covariance Ez, and then aim to recover the graph Ag by lackling the aforementioned network inverse problem. In this paper, we propose a fresh learning- solution using training examples T = ( 10
Remark (Graph convolutional data model in context). To further elaborate on the relevance and breadth of applicability of the graph convolutional (or network diffusion) signal model x H (AL; w,, we would like to elucidate connections with related work for graph structure learning. Note that while we used the diffusion-based generalive model to formulate the problem above, we do not need it as an actual mechanistic process. Indeed, like in eq. ! the only thing we ask is for the data covariance Ao = Eg to be some analytic function of the latent graph Ar. This is not extrancous to workhorse statistical methods for network topology inference, which (implicitly) make specific choices for these mappings, e.g. (i) correlation networks [Klaczyk] 2009, Ch. 7) rely on the identity mapping E, = A (ii) Gaussian graphical model selection methods, such as graphical lasso in ς oYuan & Lin 2007 Friedman al. 2008), adopt Σ, A and (ili) undirected structural equation models = ALæ + w which implies Σg [I Az)-² Mateos et al. 2019). Accordingly, these models are all subsumed by the general framework we put forth here. For recent and inspiring supervised learning approach to Gaussian graphical model selection rooted on graphical lasso, the interested reader is referred to (Shrivastava et al. 2020).
Network deconvolution and denoising. The network deconvolution problem is to identify a sparse adjacency matrix At that encodes direct dependencies, when given an adjacency matrix Ao containing extrancous indirect relationships. The problem broadens the scope of signal deconvolution to networks and can be tackled by attempting to invert the mapping Ao = Ar. (I Ar,) = Ar, + A + A ...... This solution proposed in (Feizi et al., 2013) assumes polynomial relationship as in eq. [.. but for the particular case of single-pole, singlo-zcro graph filter with very specific filter coefficients jcf. eq. (1) with ho 0 and h2 = 1, k 2 1]. This way, the indirect dependencies observed in Ao arise due to the higher-order convolutive mixture terms A + A2, +... superimposed to the direct interactions in As. we wish to recover. Here we adopt a general data-driven learning approach in assuming that Ao can be written as a polynomial in Ar, but being agnostic to the form of the filter, thus broadening the problem setting in (Feizi et al.) 2013). Unlike the problem outlined in the previous subsection, here Ao need not be a covariance matrix. Indeed, Ao could be a corrupted graph we wish to denoise, obtained via an upstream graph learning mothod. A related but different deconvolution problem was put. forth in (Li al. 2021]. Rather than recovering graph structure, the goal therein is to reconstruct the input graph signals from smoothed nodal representations generated by graph convolutional network (GCN) (Kipf & Welling 2017).
Inferring structural brain networks from functional MRI (IMRI) signals. Brain connectomes encompass networks of brain regions connected by (statistical) functional associations (FC) or by anatomical white matter fiber pathways (SC). The latter can be extracted from timo-consuming tractography algorithms applied to diffusion MRI (dMRI), which are particularly fraught due to quality issues in the data (Ych et al,, 2021). FC represents pairwise correlation structure between Mbood--oyge---eve-ddeeeeded (BOLD) signals
4

---

Published in Transactions on Machine Learning Research (06/2023)
measured by fMRI. Deciphering the relationship between SC and FC is very active area of research (Abdelnour et al. 2014: Honey et al. 2009) and also relevant in studying neurological disorders, since it is known to vary with respect to healthy subjects in pathological contexts (Gu et al., 2021). Traditional approaches exploring the SC-FC coupling go all the way from correlation studies (Greicius et al. 2008) to large-scale simulations of nonlinear cortical activity models (Honey et al., p009). More aligned with the problem addressed here, recent. studies have shown that linear diffusion dynamics can reasonably model the SC-FC coupling (Abdelnour et al. 2014). Using our notation, the findings in (Abdelnour et al., 2014) suggest that the covariance Ao = Σx of the functional signals (i.e., the FC graph) is related to the the sparse SC network AL via the model in eq. [ Similarly, (Liang & Wang 2017) contend FC can be modeled as weighted sum of powers of the SC matriz, consisting of both direct and indirect effects along varying paths. There is evidence that FC links tend to exist where there is no or little structural connection (Damoiseaux & Greicius 2009), a property naturally captured by eq. [ These considerations motivate adopting our graph strcuture learning method to infer SC patterns from iMRI signals using the training set T == (FC SS (Section 5.2 a significant problem for several reasons. The ability to collect only FC and get informative estimates of SC open the door to large-scale studies, previously constrained by the logistical, cost, and computational resources needed to acquire both modalities.
4 Graph Deconvolution Network
Here we present the proposed GDN model, a NN architecture for graph structure learning that we train in a supervised fashion. In the sequel, we obtaín 'conceptual' iterations to tackle an optimization formulation of the network inverse problem (Section 4.I), unroll these iterations to arrive at the parametric, differentiable GDN function (Ao;; ①) we train using graph data, T (Section 4.2 and describe architectural customizations to improve performance (Section 4.3
4.1 Iterative optimization as NN architectural blueprint
Going back the inverse problem of recovering a sparse adjacency matrix Ar from the mixture Ao in eq. [1 if the graph filter HAA:h) were known - but recall it is not we could attempt to solve
A || + 1 2 H(A:)))
ÅL E arg min
(3)
AEA
where A > 0 trades off sparsity for reconstruction error. The convex set A = (A E NNN (A) = 0, A, = A. 2 0, Vi.j E [V} encodes the admissibility constraints on the adjacency matrix of an undirected graph: hollow diagonal, symmetric, with non-negative edge weights. The & norm encourages sparsity in the solution, being a convex surrogate of the edge-cardinality function that counts the number of non-zero entries in A [Tibshirani 1996). Since Ao is often noisy observation or estimate of the polynomial H(A;; h), it is prudent to relax the equality in eq. ( and minimize the squared residual errors instead. As described carlier, this could be the case when Ao Ex (e.g., brain FC) and the empirical covariance matrices are estimated from finite data.
The composite cost in problem 3 is a weighted sum of a non-smooth function |||| and a continuously But our end goal here is not to solve problem 3 iteratively, recall we cannot even formulate the optimization problem because ((A:h) is unknown. To retain the essence of the problem structure and motivate a parametric model to learn approximate solutions from T, it suffices to settle with 'conceptual' proximal gradient (PG) iterations (k henceforth denote iterations, A E A)
differentiable, non-convex function g(A) := -
Alk+1 = ReLU (A[K) - -V -110)))))11))
(4)
for A: - 0, 1, where is a step-size parameter in which we have absorbed X. These iterations implement a gradient descent step on g followed by the l norm's proximal operator; see Appendix AA. for a detailed derivation of eq. (4) and (Parikh & Boyd] 2014) for more on PG algorithms. Due to the non-ncgativity constraints in A, the & norm's proximal operator takes the form of a 7-shifted ReLU on the off-diagonal entries of its matrix argument. Also, the operator sets diag( A lk + 1]) = 0.
5

---

Published in Transactions on Machine Learning Research (06/2023)
We cannot run the iterates in without knowledge of h, and we will make no attempt estimate h. Instead, our approach in the next section is to unroll and truncate these iterations (and in the process, convert. unknowns to learnable parameters ②) to arrive at the trainable GDN parametrie model (AAo; O). This way, the iterations in 4 will scrve as NN architectural blueprint to tackle the supervised problem stated in Section 2
4.2 Learning to infer graphs via algorithm unrolling
The idea of algorithm unrolling can be traced back to the seminal work of (Gregor & LeCun 2010 In the context of sparse coding, they advocated identifying iterations of PG algorithms with layers in a deep network of fixed depth that can be trained from examples using backpropagation. One can view this process as effectively truncating the iterations of an asymptotically convergent procedure, to yield a template architecture that learns to approximate solutions with substantial computational savings relative to the optimization algorithm. Beyond parsimonious signal modeling, there has been surge in popularity of unrolled deep networks for wide variely of applications; see e.g., (Monga el al. 2021 | for a recent tutorial treatment. focused on signal and image processing.
However, to the best of our knowledge this approach is yet to be fully explored for graph structure learning. L2G (Pu al. 2021) and GL.A (Shrivastava al. 2020) represent recent inspiring attempts, but neither are NNs (mcaning compositions of affine layers followed by point-wisc nonlincar activations). Both process the input data and the hidden representations in nonlinear fashion (pre activation) involving, e.g., cross- parameter products or parameter inverses within their layers, challenging optimization during training (Monga ct al. 2021 ). GLAD includes a matrix square root operation which poses scalability issues due to its cubic complexity. A detailed account of the differences between GDN and GLAD is given in Appendix A.7
Building on the algorithm unrolling paradigm and starting from the PG iterations (4, we design a non-linear, parameterized, food-forward NN architecture that can be trained to prodict the latent graph A, = (AAo: θ). We aim a NN architecture for its favorable optimization properties, yet desire to retain expressiveness. To this end, we approximate the gradient Vg(A) (A) by retaining only linear terms in A, and build deep network by composing layer-wise wise linear filters and point- boint nonlinearites to capture higher-order interactions in the generative process HA; == >x k= A4 This rationalc is similar to the one followed in developing the GCN model (Kipl & Welling 2017 Therein, it is argued that a rich class of convolutional filter functions (e.g., those parameterized via. Chebyshev polynomials) can be recovered by stacking multiple linear GCN layers.
In more detail, we start by simplifying (g(A) (derived in Appendix A.2 by dropping all higher-order terms in A, namely
k-1
Vg(A)=
hi
A - AoA" + V V Th (A: h)
23 h1 Ao - h2(AAA + AAo) + (2hg/2 + h2)
(5)
Notice that V Tr hA) is polynomial of degree 2K - 1. Hence, we keep the linear term in A but drop the constant offset that is proportional to the identity matrix I, which is inconsequential to adjacency matrix updates with null diagonal. An affine approximation leads to more benign optimization landscapes when it comes to training the resulting GDN model; see Section [ for ablation studies exploring the effects of higher-order terms in the approximation. All in all, the simplified PG iterations become
Alk 1] = ReLU (a + 3(Ao A + Alm)Ao) + + - -711
(6)
where A[0] E A and we defined a == (1 2T/p/ Th3), B = Th2, and 7 == Th1- The latter parameter triplet encapsulates filter (i.e., mixture) coefficients and the X-dependent algorithm step-size T, all of which are unknown in practice.
The GDN architecture is thus obtained by unrolling the iterations 6 into deep NN; see Figure ! This entails mapping each individual iteration into layer and stacking prescribed number D of layers together
6

---

Published in Transactions on Machine Learning Research (06/2023)
[
Single layer
[A] + 3 [K] + Ao)
A lk +1
ReLU
Ao
2
[0)
A
GDN
Layer 1
Layer 2
Layer D
Ao
Figure 1: Schematic diagram of the GDN architecture obtained via algorithm unrolling.
to form ((Ao;)) The unknown filter coefficients are treated as learnable parameters θ == {a, 3, r) which are shared across layers as in recurrent neural networks (RNNs), The reduced number of parameters relative most typical NNs is a characteristic of unrolled networks (Monga et al., 2021). A first noteworthy property of the GDN layer in eq. 6b is that it preserves symmetry, namely if the input A to layer k + 1 is equivariant to node
symmetric matrix then is the output Alk 1]. Second, GDN paramctrizations permutations essential property that guarantees the learnt graph topologies are not alfected by the order
in which we label the vertices.
Proposition 1 (Permutation equivariance of GDNs). Let Pe {0, 10 be a permutation matrix and consider identical row and column permutations of the observed graph Åo = pT Ao as well as the initial state A[0] P A Let ((Ao; O A(O]) denote the GDN output corresponding to input Ao, where for convenience make explicit that the model is parameterized by the initial state A[O]. Then, it holds that GDNs are permutation equivariant, namely
& Ao: θ, A(o) = PT θ, A(0)P.
(7)
Proof. Let Alk + 1] [[A[k), Ao) represent the layer k + mapping in eq. 6 Then by regrouping terms, recalling that P p=p=1 I, and using the fact that ReLU(.) and 11 are both permutation equivariant, we find via simple algebraic manipulations that f(PT A[k),, = P [[AAk]]AA)).. the permutation equivariance property, the GDN architecture ((Ao: A(0) that is a composition of D such layers is also permutation equivariant; meaning that eq. (7) holds true.
Because each layer satisfies
In the next section, we will explore a few customizations the vanilla GDN architecture in order broaden the model's expressive power. Given a training set T = {40) (0) A (T learning is accomplished by using mini-batch stochastic gradient descent to minimize the task-dependent loss function L(@) in eq. ( We adopt a hinge loss for link prediction and mean-squared, /absolute error for the edge-weight regression task. For link prediction, we also learn threshold & R. to binarize the estied edge weights and declare presence or absence of edges; see Appendix A.5 for all training-related details.
The iterative refinement principle of optimization algorithms carries over to our GDN model during inference. Indeed, we start with an initial estimate A() E A and use a cascade of D linear filters and point-wise nonlinearities to refine it to an output A, = & (oo θ); see Fig. for an example drawn from the experiments in Section 5.1 Matrix A is a hyperparameter we can select to incorporate prior information on the sought latent graph, or it could be learned as it is customary with the intial state in RNNSs; see Section 4.3 The input graph Ao which we aim to deconvolve is directly fed to all layers, and contributes to defining non-uniform soft thresholds 7Ao T11T which sparsify the layer's output. One also interpret aA + ((An AAo) as a first-order, symmetry- preseervvigreegngggggg filter defined on graph Ao. which is used to process A E A - here viewed as structured graph signal with N features per node to invoke this GSP insight. In its simplest rendition, the GDN leverages clements of RNNs and GCNs Kipf & Welling 2017
4.3 GDN architecture customizations
Here we outline several customizations and enhancements to the vanilla GDN architecture of the previous section, which have empirically found to improve graph structure learning performance.
7

---

Published in Transactions on Machine Learning Research (06/2023)
Ao
A[1]
A
ÅL AL
Intermediate outputs of a GDN model trained on the random geometric (RG) graphs in Table ! GDNs refine the input Ao to more closely resemble A,, as we go deeper in the network. Notice how the predicted output adjacency matrix A recovers the connectivity structure in AL.
Figure 2:
Incorporating prior information via algorithm initialization. By viewing our method as an iterative refinement of an initial graph A [0), one can think of A as a best initial guess, or prior, over At. A simple strategy to incorporate prior information about some edge (i.j), encoded in A, that we view as a random variable, would be to set A[O],} [[As]]- This technique is adopted when training on the HCP-YA dataset. in Section 5 by taking the prior A[O] to be the sample mean of all latent (i.e. SC) graphs in the training set T. This encodes our prior knowledge that there are strong similarities in the structure of the human brain across the population of healthy young adults. When AL is expected to be reasonably sparse, we can set A = 0, which is effective as we show in Table [ Recalling the connections drawn between GDNs and RNNss then the prior [ plays a similar role to the initial RNN state and thus it could be learned (Hinton 2013 In any case, the ability to seamlessly incorporate prior information to the model through the initial state an attractive feature of unrolling-based NNs (Monga et al., 2021 and of GDNs in particular.
Multi-Input Multi-Output (MIMO) filters. So far, in each layer have single learned filter, which takes an N x N matrix as imput and returns another N x N matrix at the output. After going through the shifted ReLU nonlinearity, this refined output adjacency matrix is fed as input to the next layer; a process that repeats D times. More generally, we can allow for multiple input channels (i.c., a tensor), as well as multiple channels al the output, by using the familiar convolutional neural network (CNN) methodology. This way, each output channel has its own filter parameters associated with every input channel. The j-th output channel applics its lincar filters to all input channels, aggregating the results with reduction operation (e.g., mean or sum), and applies a point-wise nonlinearity (here a shifted ReLU) to the output. This allows the GDN model to learn many different filters, thus providing richer learned representations and effectively relaxing the simplifying assumption we made regarding a common filter across graph pairs in eq. [ Notice that the aforementioned enhancement to the GDN architecture does not compromise its permutation equivariance property. Full details of MIMO filters are in Appendix A.4
Decoupling layer parameters. Thus far, we have respected the parameter sharing constraint imposed by the unrolled PG iterations. We now allow each layer to learn a decoupled MIMO filter, with its own set of parameters mapping from C.., input channels to Ch. output channels. As the notation suggests, C and Che need not be equal. By decoupling the layer structure, we allow GDNs to compose different learned filters to create more abstract features (as with CNNs or GCNs). Accordingly, it opens up the architect ural design space to broader exploration, e.g., wider layers early and skinnier layers at the end. Exploring this architectural space is beyond the scope of this paper and is left as future work. In experiments, we use GDNs for which intermediate layers E {2...., D 1} have = C Ch i.c., a flat architecture. We denote models with layer-wise shared parameters and independent. parameters as GDN-S and GDN, respectively.
5 Experiments
We present comprehensive experiments on link prediction and edge-weight. regression tasks using synthetic data (Section 5.1 as well as real HCP-YA neuroimaging and social network data (Section 5.2).
Graph generation. In the synthetic data experiments we consider three test cases whereby the latent graphs are respectively drawn from ensembles of Rényi (FR), random geometric (RG), and Barabási-A A Albert.
8

---

Published in Transactions on Machine Learning Research (06/2023)
(BA) random graph models; see e.g., Kolaczyk, 2009
We study an additional scenario where we use SCs from HCP-YA (referred to as the 'pseudo-synthetic" case because the latent graphs are real structural brain networks, while the signals used to estimate covariance matrices Ao = Σg are synthetic; see Section 5.1
Figures of merit. In the link prediction task, performance is evaluated using error := incorrectly predicted edges total possible edges
For regression, we adopt the mean-squared-crror (MSE) or mean-absolute--rror (MAE) as figures of merit.
Baselines. We compare GDNs against several relevant graph structure learning baselines: Network Deconvo- lution (ND) which uses a spectral approach to directly invert a very specific convolutive mixture (Feizi et al, 2013); Spectral Templates (SpecTemp) that advocates a convex optimization approach to recover sparse graphs from noisy estimates of Ao's eigenvectors (Segarra et al.] 2017); Graphical LASSO (GLASSO), a regularized MLE of the precision matrix for Gaussian graphical model selection (Friedman et al. 2008); a learned unrolling of altornating minimization (AM) on the GLASSO objective (GLAD) (Shrivastava et al. 2020); a learned unrolling of a primal-dual splitting (PDS) algorithm to minimize a graph recovery criterion which promotes a signal smoothness prior (I.2G) (Pu et al. )021); a GNN encoder-decoder model on graph Ao using the degree vector Aol as nodal signals, with encoders based on the GCN (GCN-D) or the graph isomorphism network (Xu ct al., 2019) (GIN-D); least-squares fitting of h followed by non-convex optimization to solve problem 33 (LSOpt); and Hard Thresholding (Threshold) to assess how well simple cutoff rule can perform. GLASSO, ND, and SpeeTemp generally tune their regularization parameters in an unsupervised way, but here they are tuned with supervision; see A.. for details on all baselines we implemented. The network inverse problem only assumes knowledge of the observed graph Ao for inference - there is no nodal information available and GDN layers do not operate on node features; see Appendix [A.4] for a discussion on how GDNs can be extended to handle node fcatures. As such, popular GNN-based methods used in the fink-prediction task (Wang et. al. 2019 Kazi el al. 2023 Veličković et al.] 2020 Kipl et al.) 2018 Kipl & Welling 2016 Zhang & Chen (018)) are not directly applicable here. We carefully reviewed the network inference literature and chose the most performant basclincs available for this problem sctting.
GDN-related hyperparameters and experimental setting. Unless otherwise stated, in all the results that follow we use GDN(-S) models with D = & layers, C 8 channels per layer, take prior A[0] = 0 on all domains except the SCs where we use the sample mcan of all SCs in the training set, and train using the ADAM optimizer (Kingma & Ba 2015) with learning rate of 0.01 and batch size of 200. We use one Nvidia T4 GPU; models take < 15 minutes to train on all datasets. detailed account on the GDN hyperparameter search process is included in the supplementary material; see Appendix A.8.
5.1 Learning graph structure from diffused signals
A set of latent graphs are either sampled from RG, ER, or the BA model, or laken as the SCs from the HCP-YA dataset. In an attempt to make the synthetic latent graphs somewhat comparable to the 68 node
SCs, we sample connected N = 68 node graphs with edge densities in the range [0.5, 0.6] when fcasible- -
(i) typical values observed SC graphs. To generate each observation A we simulated P 50 standard
Normal white signals diffused over A from which we form the sample covariance E, as described in Section The value of P is chosen based on study of robustness to noise in the observations across multiple models, which showed only marginally improved performance when averaging more signals as shown in Table 2 With respect to the diffusion graph filters we let K = 2, and sample the filter coefficients h E R in (A h) uniformly from the unit sphere. To examine robustness to the realizations of h, we repeat this data generation process three times (resampling the filter coefficients). We thus create three different datasets for each graph domain (12 in total). For the sizes of the training/ validation//ee /eestoot/ttttttt splits, the pseudo-synthelic domain uses 913 /50/500/100 and the synthetic domains use 113/500/500/
Table 1 tabulates the results for synthetic and pscudo-synthetic experiments. For graph models that exhibit localized connectivity patterns (RG and SC), GDNs significantly outperform the baselines on both lasks. For the SC test case, GDN (GDN-S) reduces error relative to the mean prior by 27.5 + 1.7% (23.0 + 1.7%) and MSE by 37.3 = 0.8% (23.2 + 0.5%). Both GDN architectures show the ability to learn such local patterns, with the extra representational power of GDNs (over GDN-S) providing an additional boost in performance. All models struggle on BA and ER with GDNs showing better performance even for these cases.
9

---

Published in Transactions on Machine Learning Research (06/2023)
Table 1: Mean and standard error of the test performance across both tasks (Top: link-prediction, Bottom: odge-weight regression) on each graph domain. Bold denotes best performance.
Models
RG
ER
BA
SC
GDN
41.9
27.5 2te-
8.9 :2-2
GDN-S
5.51
40.8 40.8 Sat
27. 97.64
9.4
L2G
44.3
43. 9
34.9+1ca
31.7+ 3
GCN-D
30. 30.723e-4
41.8
27.6
31.6 11.62te-3
Error (%)
GIN-D
20.0+20-8
43. 23.68
28.3 9o-4
33. 33.3++.. 8c-8
GLAD
6.3
43.7 13.723mm 33e11
35-02se-2
11.81 2e-3
GLASSO
8.8
43.2 13.2
34.9 9e-8
20.0
ND
9.4+3e-1
43.9 43.9+1e2
34.1+ 8e-8
21.3+
SpecTTemy
11.1 3e-1
44.4
30.21 2e-1
30.011e1
LSOpt
24.2255-0
42.5
28.0
31.51se-3
Threshold
26.8
42.9
32.3
21.7
GDN
4.26-2 2 +4
2.3e-1+20
1.8c- 1 +2c-3
5.3e-3
GDN-S
6.0e-2+2e4
2.3e 1
2.7e-1 22-22
6.5e-3+
L2G
2.5c-1 4
2.5c- 1
2.3c- +
4.9c-2
GCN-D
2. 1
2.4e- 1 136-
2.0e- 1 1ate
1.1e- 1
MSE
GIN-D
2.0e- 1
2.6e- 1+c..8
1.9e- 1.9e-1+7e4. 1+7444
2. 2.6e- 1+5e.a
GLAD
7.8c 7 #te-
2.4c- 10-114402 44-3
1.9c- + 4e-3
1.40-2
GLASSO
2.0e- 1 ++30.3
2.8e-1+2e2
2.6e- +2:22
ND
1.8e- 1 22003
2.4e- 155-1
2.20- 1 11:-3
1.60222Te
SpecTemm LSOpt
5. 2
5.3e- 1
3.3e-
2
1.5e- 1 4
9.9c-2+2e.
2.5e- 22203
2.0c- 1 1
6.1e-0
Robustness to noise. To evaluate the robustness of GDN's to noise the observations Ao, we train a GDN model with D = 30 and C' on multiple datasets with varying levels of noise. We impliclty control the level of noise by changing the number of signals P used to estimate Ao taken to be the sample covariance matrix Σ.. The data are the same in Table except here we perform diffusions using one of the sampled filter coefficients h E R The results of this experiment on both GDN and Threshold are reported in Table 2 We find that GDN's are robust to noisy observed graphs Ao, achieving less than 10% error with only 20 observed signals. Up to P= 50, GDN performs significantly better when given more signals, with only marginal increases in performance after this point. This experiment also highlights the advantage of MIMO filters as a non-MIMO GDN of depth 30 is outperformed by the MTMO GDN of depth on the P= case.
Scaling and size genoralization: Deploying larger graphs. GDNs learn the parameters of graph convolutions for the processing of graphs making them inductive: we deploy the learnt model on larger graph size domains. Such a deployment is feasible on larger graphs due to the O(N2) time and memory complexity of GDNs; the same is true for L2G. This stands in contrast to the O(N3) time complexity of SpecTemp, GLASSO, ND, and GLAD. SpecTemp and ND perform a full eigendecomposition of Ao, the iterative algorithms for solving GLASSO incur cubic worst-case complexity (Zhang et al. 2018), and GLAD performs matrix square root each layer. GLAD also requires deeper unrolled networks, discounted intermediate losses, and multiple embedded MLPs, which have significant practical implications on memory usage; also Appendix A.7 To substantiate this qualitative discussion, in Appendix A.& we have included supplementary experimental runtime and memory usage analyses showing GDN compares favorably to SpecTeem, L2G, and GLAD.
To test the extent to which GDNs generalize when N grows, we trained GDN(-S) on RG graphs with size = 68, and tested them on RG graphs of size N [68, 75, 100, 200, 500, 1000, 2000), with 200 graphs of each size. As graph sizes increase, we require more samples in the estimation of the sample covariance to maintain constant signal-to-noise ratio. To simplify the experiment. and its interpretation, disregard estimation
10

---

Published in Transactions on Machine Learning Research (06/2023)
Table 2: Mean and standard error of the test performance on link prediction task for a varied number P of sampled signals N = 68 RG graphs.
Models
5
20
50
100
200
1000
GDN
17.9
9.6
6.818m4
6.418m4
6.3
Threshold
31.82 28
27.1
26.8126-3
26.817m-1
26.817m-1
26.817m-1
and directly use the ensemble covariance Ao = En as observation. As before, we take training/validation /allidaiionnnn split of 913/500. Figure 3 shows GDNs effectively generalize to graphs orders of magnitude larger than they were trained on, giving up only modest levels of performance as size increases. Note that the best performing baseline in Table D trained and tested on the original N = 68 domain is not comparable with GDN-S in terms of performance until GDN-S is lested on graphs almost an order of magnitude larger than those it was trained on. The GDN-S model showed a more graceful performance degradation, suggesting that GDNs without parameter sharing may be using their extra representational power to pick up on finite-size effects, which may disappear as N increases. The shared parameter constraint acts as regularization, we avoid over-fitting a given size domain to better generalize to larger graphs.
A blation studies. The choice of prior can influence model performance, as well as reduce training time and the number of parameters needed. When run on stochastic block model (SBM) graphs with N = 21 nodes and equally-sized communities (within block connection probability of 0.6, and 0.1 across blocks), for the link prediction task GDNs attain an error of 16.8 2.7e-2% 16.0 2.1.2%, 14.5 + 1.0e-2%, 14.3 = 8.8 using zeros, ones, block diagonal (using the SBM edge connectivity probabilities). and learned prior, respectively. The performance improves when GDNs are given an informative prior (here block diagonal matrix matching the graph communities) with further gains when GDNs are allowed to learn A[O].
We also study the effect of gradient truncation. To derive GDNs we approximate the gradient (g(A) by dropping all higher-order terms in A (K- = 1). The case of K = 0 corresponds to further dropping the terms lincar in A, leading to PG itcrations Alk + 1] [ALU((([]] + yAo - -711 [cf. eq. (6). We run this simplified iterative soft-thresholding model with D 8 layers and C = 8 channels per layer on the same RG graphs in Table 1 Lacking the linear term that facilitates information aggregation in the graph, the model is not expressive enough and yields higher error (MSE) of 25.7 + 1.3e-2% (1.7e-1 + 4.7e-4) for the link-prediction (edge weight regression) task. Models with K 2 2 result in unstable training with highly fluctuating loss values, which we attribute to a challenging optimization landscape as polynomial minimization is known to be notoriously hard, leading to less well-behaved gradients. This motivates our choice of K = in GDNs, which are nonetheless flexible enough to capture a rich class of (higher-order polynomial) convolutional processes by stacking multiple linear layers as in Figure]
5.2 Real Data
HCP-YA neuroimaging dataset. HCP represents a unilying paradigm to acquire high quality neuroimag- ing data across studies that enabled unprecedented quality of macro-scale human brain connectomes for analyses in different domains (Glasser et al. 2016). We use the dMRI and resting state fMRI data from the HCP-YA dataset (Van Essen et al. 2013 consisting of 1200 healthy, young adults (ages: 22-36 years). The SC and FC are projected on a common brain atlas, which is a grouping of cortical structures in the brain to distinct regions. We interpret these regions as nodes in a brain graph. For our experiments, we use the standard Desikan-Killiang atlas (Desikan et al.) 2006) with N = 68 cortical brain regions. The SC-FC coupling on this dataset is known to be the strongest in the occipital lobe and vary with age, sex and cognitive health in other subnetworks (Gu et al.) (021).. Under the considcration of the variability in SC-FC coupling across the brain regions, we further group the cortical regions into 4 neurologically defined 'lobes': frontal, parietal, temporal, and occipital (the locations of these lobes in the brain are included in Fig. 0 in Appendix A.6). We predict SC, derived from dMRI, using FC, constructed using BOLD signals acquired with fMRI. In the literature, the forward (SC to FC) problem has been mostly attempted, a non con-trivial yet simpler task due to the known reliance of function on structure in many cortical regions. We tackle the markedly harder
11

---

Published in Transactions on Machine Learning Research (06/2023)
8
7
CITO (%)
6-
GDN
4-
GDN-S
68
100
200 N (log)
500
1000
2000
4
6
Figure 3: Size Generalization: GDNs maintain perfor-
Figure Reduction in MAE (%) for different lobes; largest improvements oceur in lemporal and frontal lobes.
mance on RG graphs orders of magnitude larger than the N 68 training set.
inverse (dcconvolution) problem of recovering SC from FC, which has so far received less attention and can have major impact as described in Section Notably, this network deconvolution problem represents an ideal fit to the polynomial model 1 relating FC with SC (Liang & Wang, 2017: Abdelnour et al. 2014), plus the supervised setting is well motivated given the availability of functional and structural connectomes from the HCP-YA repository.
From this data, we extracted a dataset of 1063 FC-SC pairs, T (FC (CC)) 109 and use train- nug/validaaionn/estt split of 913/50/ 100. Taking the prior A[0] as the edgewise mean over all SCs in the training split Terain: A[0];; = mean {SC ....,,.......................... (c.... and using it directly as a predictor on the test set (tuning a threshold with validation data), we achieve strong performance on ink-prrdiction and edge- weight regression lasks on the whole brain (error 12.3%, MAE 0.0615). At the lobe level, the prior achieves relatively higher accuracy in occipital (error 1.3%, MAE 0.064) and parietal (error 6.6%, MAE 0.07) lobes compared to temporal (error 11.0%, MAE 0.053) and frontal (error = 10.5%, MAE 0.062) lobes; behavior which is unsurprising as SC in temporal and frontal lobes are affected by ageing and gender related variability in the dataset (Zimmermann et al. 2016). GDNs reduced MAF by 7.6%. 7.1%, 1.6%, and 1.3% in the temporal, frontal, parietal, and occipital networks respectively and 8.0% over the entire brain network, all relative to the already strong mean prior. The four lobe reductions are visualized in Figure We observed the most significant gains over temporal and frontal lobes; clearly there was smaller room improve performance over the occipital and frontal lobes. The resulting MAEs for the entire brain network for L2G, GLAD, GCN-D, and GIN-D are 2.23e-1, 2.07e-1, 2.23e-1, 2.10e-1 respectively; none were able to reduce MAE (or error).
In summary, our 'psoudo-synthetic' experiments in Section 5.1 show that SCs are amenable to learning with GDNs when the SC-FC relationship satisfies eq. 11 a reasonable model given the findings (Abdelnour et al.) 2014 Liang Wang, 2017). In general, SC-FC coupling can vary widely across both the population and within the sub-regions of an individuals brain for healthy subjects and in pathological contexts. When trained the HCP-YA dataset, GDNs exhibit robust performance over such regions with high variability in SC-FC coupling. Therefore, our results HCP-YA dataset justify GDNs as the leading baseline model for comparisons in this significant FC to SC task. It could also potentially serve as a baseline for characterizing healthy subjects in pathology studies in future work, where pathology could be characterized by specific deviations from our results.
Friendship recommendation from physical co-location networks. Here we use GDNs to predict Facebook ties given human co-location data. GDNs are well suited for this deconvolution problem since one view friendships encounters in addition to face-to-face contacts with friends. In terms of impact, trained model could then be useful to augment friendship recommendation engines given co-location (behavioral) data. Granted, this would also pose interesting privacy questions that are beyond the scope of this paper, and certainly call for tighter control co-location information.
direct ties, whereas co-location edges include indirect relationships due to casual
The Thiers13 dataset Géénos & Barrat 2018) monitored high school students, recording: (1) physical co-location interactions with wearable sensors over days; and (ii) social network information via survey.
12

---

Published in Transactions on Machine Learning Research (06/2023)
From this we construct dataset T of graph pairs, each with N 120 nodes, where Ao are co-location networks, i.e, weighted graphs where the weight of edge (ī,j) represents the number of times student i and j came into physical proximity, and Ar is the Facebook subgraph between the same students; additional data pre-processing details are in Appendix A.6 We trained a GDN of depth D = 11 without MIMO filters (C = 1), with learned A using a training/ validation /est split of 5000//000/11000
We achieved a test. error of 8.9 1.5e-%% a 13.0% reduction over next best performing bascline, suggesting we learn a latent. mechanism that the baselines (even GLAD, a supervised graph learning method) cannot. The performance of the baselines are follows: Threshold (10.2 + 1.2e-2%), ND (10.2 1.2e-2%).. GLASSO (NA could not converge for any of the a values tested), SpecTemp (10.7 + 6.9e-2%), GLAD (14.7 + 1.3e-2%) GCN-C (11.7 + 3. 3.7,,),), and GIN-D (11.5 + 4.0e-4%). See Appendix A.6 for further details.
6 Conclusions
In this work we proposed the GDN, an inductive model capable of recovering latent graph structure from observations of its convolutional mixtures. Bringing bear the principle of algorithm unrolling, we introduce novel deep learning solution to network inverse (deconvolution) problem in supervised fashion, that is well motivated e.g. to predict whole brain structural connectivity from functional signals. By minimizing task-dependent loss function, GDNs learn filters to refine initial estimates of the sought latent graphs layer by layer. The unrolled NN can seamlessly integrate domain-specific prior information about the unknown graph distribution. GDNs enjoy a node permutation equivariance property, which is key to ensure the learnt graph topologies are not affected by arbitrary vertex labelings. Moreover, because GDNs: (i) are differentiable functions with respect their paramcters as well as their graph input; and (ii) offer explicit control on complexity (leading to fast inference times); one can envision GDNs as valuable components in larger (even online) end-to-end graph representation learning systems. This way, while our focus here has been exclusively for graph structure learning, the impact of GDNs can permeate broader graph inference tasks.
References
Farras Abdelnour, Henning U. Voss, and Ashish Raj. Network diffusion accurately models the relationship between structural and functional brain connectivity networks. Neuroimage, 90:335-333 347, Apr. 2014.
Anaconda. Anaconda Software Distribution, 2020. URL https://doca anaconda.
Alain Barrat, Mare Barthelemy, and Alessandr Vespignani. Dynamical Processes on Compler Networks. Cambridge University Press, New York, NY, 2008.
Michael M Bronstein, Joan Bruna, Yann LeCun, Arthur Szlam, and Pierre Vandergheynst. Geometric deep lcarning: Going beyond Euclidean data. IEEE Signal Process. Mag., 34(4):18-42, -2, 2017.
Jessica S Damoiseaux and Michael D Greicius. Greater than the sum of its parts: A review of studies combining structural connectivity and resting-state functional connectivity. Brain Struct. Func., 213(6): 525-533, 2009.
Hadi Daneshmand, Manuel Gomez-Rodriguez, Le Song, and Bernhard Schoelkopf. Estimating diffusson network structures: Recovery conditions, sample complexity soft-thresoldding algorithm. Tn International Conference on Machine Learning, pp. 793 801. PMLR, 2014.
Morris H. DeGroot. Reaching a consensus. J Am Stat Assoc., 69:118-121, 1974.
Arthur P Dempster. Covariance selection. Biometrics, pp. 157 175, 1972.
Rahul S Desikan, Florent Ségonne, Bruce Fischl, Brian T Quinn, Bradford Dickerson, Deborah Blacker, Randy I, Buckner, Anders M Dale, Paul Maguire, Bradley Hyman, et al. An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest. Neuroimage, 31 (3):968 1-968-900, 2006.
13

---

Published in Transactions on Machine Learning Research (06/2023)
Xiaowen Dong, Dorina Thanou, Pascal Frossard, and Pierre Vandergheynst. Learning Laplacian matrix in smooth graph signal representations. IEEE Trans. Signal Process., 64(23):6160-61)) 6100,,17, 2016.
Xiaowen Dong, Dorina Thanou, Michael Rabbat, and Pascal Frossard. Learning graphs from data: A signal representation perspective. IEEE Signal Process. Mag., (3)3):44 63, 2019.
Mathias Drton and Marloes H Maathuis. Structure learning graphical modeling, Annu. Rev. Stot. Appl, 4: 365 393, 2017.
Vijay Prakash Dwivedi, Anh Tuan Luu, Thomas Laurent, Yoshua Bengio, and Xavier Bresson. Graph neural networks with learnable structural and positional representations. In International Conference on Learning Representations, 2022. URL http::///ggpeeeeeeee.e.ee//rrr/a
Soheil Feizi, Daniel Marbach, Muriel Medard, and Manolis Kellis. Network deconvolution as a general method to distinguish direct dependencies in networks. Nat. Biotechmol, (1())7726 733, 2013.
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, ((3):432-44 441, 2008.
Fernando Gama, Joan Bruna, and Alejandro Ribeiro. Stability properties of graph neural networks. IEEE Trans. Signol Process., 68:5680 5695, 2020a.
Fernando Gama, Elvin Isufi, Geert Leus, and Alejandro Ribeiro. Graphs, convolutions, and neural networks: From graph filters to graph neural networks. IEEE Signal Process. Mug., 37(6):128-11) 138, 2020b.
Mathieu Génois and Alain Barrat. Can co-location be used as proxy for face-to-face contacts? EPJ Dala Science, 7(1):11, May 2018.
Georgios B. Giannakis, Yanning Shen, and Georgios Vasileios Karanikolas. Topology identification and learning over graphs: Accounting for nonlinearities and dynamics. Proceedings of the IEEE, 106(5):787-80) 2018.
Matthew F Glasser, Stephen M Smith, Daniel Marcus, Jesper LR Andersson, Edward Auerbach, Timothy EJ Behrens, Timothy S Coalson, Michael P Harms, Mark Jenkinson, Steen Moeller, et al. The human connectome project's neuroimaging approach. Nat. Neurosci., 19(9):1175-1187, 2016.
Karol Gregor and Yann LeCun. Learning fast approximations of sparse coding. In International Conference on Machine Learning, pp. 399-406, 2010.
Michael Greicius, Kaustubh Supekar, Vinod Menon, and Robert Dougherty. Resting-state functional connectivity reflects structural connectivity in the default mode network. Cereb Corlez, 19:72 8, 12 2008
Zijin Gu, Keith Wakefield Jamison, Mert Rory Sabuncu, and Amy Kuceyeski. Heritability and interindividual variability of regional structure-fuuction coupling. Not. Commun, 12(1): 2021.
Geoffrey E. Hinton. Advanced machine learning: Recurrent neural networks, 2013. URL https://wwwy cs pdf
toronto.
C. Honey, O. Sporns, L. Cammoun, X. Gigandet, J. Thiran, R. Meuli, and P. Hagmann. Predicting human resting-state functional connectivity from structural connectivity. Proc. Notl. Acad. Sei. U.S.A, 106: 2035-2040, 2009.
Vassilis Kalofolias. How to learn graph from smooth signals. In Artificial Intelligence and Statistics, pp- 920 929, 2016.
Anees Kazi, Luca Cosmo, Seyed-Ahmad Ahmadi, Nassir Navab, and Michael M. Bronstein. Differentiable graph module (dgm) for graph convolutional networks. IEEE Trans. Pattern Anol. Mach. Intell., 45(2): 1606-1617, 2023.
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In International Conference Learning Representations, 2015.
14

---

Published in Transactions on Machine Learning Research (06/2023)
Thomas N. Kipf and Max Welling. Variational graph auto-encoders. CoRR, abs/1611.07308, 2016. URL
http://arxiv httpaaaxivooror/aaaa aatpabb.../..00/0
Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In
International Conference on Learning Representations, 2017.
Thomas N. Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In International Conference on Machine Learning, volume 80, pp. 2688-2697, 10-15 Jul 2018.
Eric D. Kolaczyk. Statistical Analysis of Network Data: Methods and Models. Springer, New York, NY, 2009.
Jia Li, Jiajin Li, Yang Liu, Yueting Li Jianwei Yu, and Hong Cheng. Deconvolutional networks on graph data. In Advances in Neural Information Processing Systems, 2021.
Hualon Liang and Hongbin Wang. Structure-funnttion network mapping and its assessment via persistent homology. PLOS Compul. Biol, 13(e1005325):1 19, Jan. 2017.
Gonzalo Mateos, Santiago Segarra, Antonio G. Marques, and Alejandro Ribeiro. Connecting the dots: Identifying network structure via graph signal processing, IEEE Signal Process. Mag., 36(3):16 43, May 2019.
Vishal Monga, Yuelong Li, and Yonina C. Eldar. Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing. IEEE Signal Process. Mag., 38(2):18-44, 2021.
Antonio Ortega, Pascal Frossard, Jelena Kovačević, José M. F. Moura, and Pierre Vandergheynst, Graph signal processing: Overview, challenges, and applications. Proc. IEEE, 106(5):808 828, 2018.
Neal Parikh and Stephen Boyd. Proximal algorithms. Foundations and Trends in Optimization, 133)1277-33,
2014.
Bastien Pasdeloup, Vincent Gripon, Grégoire Mercier, Dominique Pastor, and Michael G. Rabbat. Charac- terization and inference of graph diffusion processes from observations of stationary signals. IEEE Trans. Signal Inf. Process. Netw., (33):481496,, 2018.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performanc deep learning library. In H. Wallach, H. Larochelle, A. Beygelzimer, d' A'AchééBucc, E. Fox, and R. Garnett (cds.), Advances in Neural Informa- tion Processing Systems 32, pp. 8024 8035. Curran Associates, Inc., 2019. URL httpp:/papprs. http://ppaprr....ippppppp neurips .df
Xingyue Pu, Tianyue Cao, Xiaoyun Zhang, Xiaowen Dong, and Siheng Chen. Learning to learn graph topologies. In Advances in Neural Information Processing Systems, 2021.
Santiago Segarra, Antonio G. Marques, Gonzalo Mateos, and Alejandro Ribeiro. Network topology inference from spectral templates. IEEE Trans. Signal Inf. Process. Netw., 3(3):467-483, August 2017.
Harsh Shrivastava, Xinshi Chen, Binghong Chen, Guanghui Lan, Srinvas Aluru, and Le Song. GLAD: Learning sparse graph recovery. In International Conference on Learning Representations, 2020.
Daniel A. Spielman and Nikhil Srivastava. Graph sparsification by effective resistances. STAM J. Comput., 10(6):1913-1925 December 2011.
Pablo Sprechmann, Alex M. Bronstein, and Guillermo Sapiro. Learning efficient sparse and low rank models. IEEE Trans. Pattern Anal. Mach. Intell., 179)118211 2015.
Robert Tibshirani. Regression shrinkage and selection via the lasso. J R Stat Soc Series B Stat Methodol,, 58 (1):267-288) 1996.
15

---

Published in Transactions on Machine Learning Research (06/2023)
David C Van Essen, Stephen M Smith, Deanna M Barch, Timothy EJ Behrens, Essa Yacoub, Kamil Ugurbil, Wu-Minn HCP Consortium, et al. The WU-Minn Human Connectome Project: An overview. Neuroinage, 80:62 79, 2013.
Petar Veličkovié, Lars Buesing, Matthew C. Overlan, Razvan Pascanu, Oriol Vinyals, and Charles Blundell. Pointer graph networks. In Advances in Neural Information Processing Systems, 2020.
Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E. Sarma, Michael M. Bronstein, and Justin M. Solomon. Dynamic graph CNN for learning on point clouds. ACM Trans. Graph., 38(5), October 2019.
Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka, How powerful are graph neural networks? In International Conference on Learning Representations, 2019.
Chun-Hung Yeh, Derek K Jones, Xiaoyun Liang, Maxime Descoteaux, and Alan Comnelly, Mapping structural connectivity using diffusion MRI: Challenges and opportunities. J. Magn. Reson., 53(())116661666 2021.
Gilad Yehudai, Ethan Fetaya, Eli Meirom, Gal Chechik, and Haggai Maron. From local structures to size generalization in graph neural worworks. In International Conference on Machine Learning, pp. 11971- 11986. PMLR, 2021.
Ming Yuan and Yi Lin. Model selection and estimation in the Gaussian graphical model. Biometrika, 94(1): 19 35, 2007.
Muhan Zhang and Yixin Chen. Link prediction based on graph neural networks. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Information Processing Systems, volume 31. Curran Associates, Inc., 2018.
Richard Zhang, Salar Fattahi, and Somaych Sojoudi. Large-scale sparse inverse covariance estimation via thresholding and max-det matrix completion. In International Conference on Machine Learning, pp. 5766-5775. PMLR, 2018.
Tong Zhao, Yozen Liu, Leonardo Neves, Oliver Woodford, Meng Jiang, and Neil Shah. Data augmentation for graph neural networks. In AAAI Conference on Artificial Intelligence, number 12, pp. 11015-11023, 2021.
Joelle Zimmermann, Petra Ritter, Kelly Shen, Simon Rothmeier, Michael Schirner, and Anthony R McIntosh. Structural architecture supports functional organization in the human aging brain at a regionwisc and network level. Hum. Brain Mapp, 37(7):2645 2661, 2016.
A Appendix
A.1 Model Identifiability
Without any constraints on h and AL, the problem of recovering AL from Ao H( as in eq. ① is clearly non-identifiable. Indeed, if the desired solution is Ar. (with associated polynomial coefficients h), there is always at least another solution Ao corresponding to the identity polynomial mapping. This is why adding structural constraints like sparsity on AG, will aid model identifiability, especially when devoid of training examples.
It worth mentioning that eq. ① implics the eigenvectors of A and Ao coincide. So the cigenvectors of the sought latent graph are given once observe Ao, what left to determine are the eigenvalues. We have in practice observed that for several families of sparse, weighted graphs, the eigenvector information along with the constraint A, A are sufficient to uniquely specify the graph. Interestingly, this implics that many random weighted graphs can be uniquely determined from their eigenvectors. This strong uniqueness result does not render our problem vacuous, since seldomly in practice one gets to observe Ao (and hence its eigenvectors) error free.
If one were formally study identifiability of eq. 11 (say under some structural assumptions on Ag and/or the polynomial mapping), then one has to recognize the problem suffers from inherent scaling ambiguity.
16

---

Published in Transactions on Machine Learning Research (06/2023)
Indeed, if given Ao = (AL;:)) which means the pair AL and h [hg, hi,..., hx]" T is solution, then for any positive scalar a one has that a,A, and [ho, hu ha./a.....u hxx (a ]] )] is another solution. Accordingly, uniqueness claims can only be meaningful modulo this unavoidable scaling ambiguity, But this ambiguity is lifted once we tackle the problem in a supervised learning fashion our approach in this paper. The
training samples in T := (A0 )T fix the scaling, and accordingly the GDN can learn the mechanism or mapping of interest Ao + AL. Hence, an attractive feature of the GDN approach is that by using data, some of the inherent ambiguities in eq. 1 are naturally overcome. In particular, the SpecTemp approach in (Segarra et al. 2017) relies on convex optimization and suffers from this scaling ambiguity, so it requires an extra (rather arbitrary) constraint fix the scale. The network doconvolution approach (ND) in (Feizi et al. 2013) relies on lixed, known polynomial mapping, and while does not suffer from these ambiguities it is limited in the graph convolutions it can model.
All in all, the inverse problem associated to eq. (1) is just our starting point to motivate a trainable parametrized architecture A = & (Ao: θ) that introduces valuable inductive biases to generate graph predictions. The problem we end up solving is different (recall the formal statement in Section 2 because we rely on supervision using graph examples, thus rendering many of these challenging uniqueness questions less
relevant.
A.2 Technical Details of the Proximal Gradient Algorithm
Derivation of the proximal gradient iterations
4
Here we derive the proximal gradient (PG) With a > and A a convex set, introduce the
algorithm for the composite optimization problem proximal operator of a function asp(-) RNNN + evaluated at matrix M RNXN as
1 MI
((M) pro (M) arg min ((X)+
|X
(8)
XEA
2a
With these definitions, the PG updates with fixed step-size I to solve the network inverse problem given by (k= = 0, denote iterations)
3 are
A 1] = prox (A VV((()))
(9)
Notice that r should be appropriately chosen as a function of the regularization parameter À in 3
Evaluating the proximal operator efficiently is key to the success of PG methods. For our specific case of sparse graph learning where ((A) = ||Al] and the optimization variable is constrained as A A, the proximal operator Z in has entries given by
0. Zr((Ma)= = max(0, Mij - a),, otherwise.
i=
(10)
The resulting entry-wise separable nonlinear map nulls the diagonal entries of A(k: + 1], and applies non-negative noft-thresholdin operator to update the remaining entries. Recognizing the a-shifted rectified
linear unit ReLU(M,; - a) max(0, M,y a) and substituting the proximal operator expression (10) in
one arrives at the PG iterations 4 We recall that as per 110 the matrix-valued ReLU operator in sets
diag( A + 1)) = 0. Tʃ A E A, then symmetry of all subsequent iterates A + 1), h: = 0, 2, will be preserved.
Gradient calculation. Here we give mathematical details in the calculation of the gradient Vg(A) of the component. function (AA) Aoo H(A: h) Let A, Ao be symmetric N x N matrices and recall the
17

---

Published in Transactions on Machine Learning Research (06/2023)
graph filter (AA) Ex=o A² (we drop the dependency in h to simply the notation). Then
V V42 Ao H(A] VA V Tr (A AAH(A) H(A)Ao (A))
VA [Ao (T
--1
k-r-1 AoA' 1 2 Tr Ty
116
k=
k-1
A 1 AoA + (A(
k=1
r=0
= - Al - (2 + AoA) + ( + AAoA+
+ +...] ...]+ AAoA+AoA)) Hi (A),
where in arriving at the second equality we relied on the cyclic property of the trace, and H (A) is a matrix polynomial of order 2K-
Note that in the context of the GDN model, powers of A will lead to complex optimization landscapes, and thus unstable training. We thus opt to drop the higher-order terms and work with first-order approximation of Vg, namely
= ~
A2(AAA + AAo) + (2hoh2 + +1))
A.3 Incorporating Prior Information
In Section 4.3 we introduce the concept of using prior information in the training of GDNs. We do so by encoding information we may have about the unknown latent graph Ar into A[0), the initial state which GDNs iteratively refine. the Ay's are repeated instances of a graph with fixed nodes, as is the case with the SCs with respect the 68 fixed brain regions, simple strategy to incorporate prior information about some edge ALL,, now viewed as a random variable, would be A + B[A But there more that be done. We also estimate the variance Var(A Var(Ar,,)) and use it during the training of GDN, for example taking A[kk, < (EEAA,,,,) Var(An.))), or even simpler, using resampling technique and taking [ to be a random sample in the training set. By doing so, we force the GDN to take into account the distribution and uncertainty in the data, possibly leading to richer learned representations and better performance. It also would act as a form of regularization, not allowing the model to converge on the naive solution outputting the simple expectation prior, a likely local minimum in training space. We also can learn the prior A by treating all entries as parameters directly; note due undirectedness and lack of self-loops this would require N(N 1)/2 additional parameters.
A.4 MIMO Model Architecture
MIMO filters. Formally, the more expressive GDN architecture with MIMO (Multi-Input Multi-Output) filters constructed as follows. At layer k of the NN, we take three-way tensor A E RCXNXN and produce Ax+1 E C C where C the common number of input and output channels. The assumption of having a common number of input and output channels can be relaxed, as we argue below. By defining multiplication between tensors B E xC as batched matrix multiplication: [TB], jas == T,,, and tensor-vector addition and multiplication T == T [1 and [T J.. 2y T respectively for RC, all operations extend naturally.
Using these definitions, the j-th output slice of layer k is
[A = ReLUUaLaaaaee A + (AoA: + AxAo) + TY, - -,,
(11)
where -- represents the mean reduction over the filtered input channels and the parameters are a, B, CCC RÇ. We take θ = (.,8,7, for total of C (3C + 1) trainable parameters.
E
18

---

Published in Transactions on Machine Learning Research (06/2023)
Layer k
1
C
1 ... j
C
MIMO filter
A
A
IA
[A
K+T
Figure 5: MIMO Filter: Layer takes a tensor A RC X and outputs tensor A&+1 x The i-th slice (Aplas called the i-th input channel and [As+1las called the j-th output channel.
We typically have single prior matrix A(0) and are interested in predicting a single adjacency matrix ÅL. Accordingly, we construct new tensor prior A == .., the first output channel
A(0]] CxNNN and (arbitrarily) designate
prediction ÅL Φ (A θ
(Ax+1l1....
We can also allow each layer to learn a decoupled MIMO filter, with its own set of parameters mapping from Ch imput channels to Oth output channels. As the notation suggests, Oh,, and C need not be equal. Layer k now has its own set of parameters * = (at , , where a 7/ xh.. and T² e L for a total of Che x (3Cf, + 1) trainable parameters. The tensor operations mapping inputs to outputs remains basically unchanged with respect to T1 except that the filter coefficients will depend on
Processing at generic layer k is depicted in Figure 5 Output channel j will use o: B V. E R. and o E R2, to filter all imput channels i E {1, ... J. which are collected in the input tensor A E RC,, X This produces a tensor stacked filtered input channels RR XNXN After setting the diagonal elements of all matrix slices in this tensor to O, then perform a mean reduction edgewise (over the first. mode/ /dimensii))) of this tensor, producing single N N matrix. We then apply two pointwise/elemennttee this matrix: (i) subtract T (this would be the 'bias' term in CNNs): and (fi) apply point-wisc nonlincarity (ReLU). This produces an N x N activation stored in the j-th output channel. Doing so for all oulput. channels f1,.... CCh.r, produces tensor A E R Co XNNN
operations on
Layer output normalization and zeroing diagonal. The steps not shown in the main text are the normalization steps, practical issue, and the setting of the diagonal elements to be O, projection onto the A of allowable adjacency matrices. Define U J.... = = + (AoAA + A,Ao) + + Ao E N the j-th slice in the intermediate tensor U R. used in the filtering of Ag- Normalization is performed by dividing each matrix slice of U by the maximum magnitude element in that respective slice: U..../ max U....-
Multiple normalization metrics were tried on the denominator, including the 99th percentile of all values in U U.... the Frobenius norm [U.,..le, among others. None seemed to work as well as the maximum magnitude clement, which has the additional advantage of guarantecing entries to be in [0, 1] (after ReLU), which matches nicely with: (i) adjacency matrices of unweighted graphs; and (ii) makes it easy normalize edge weights of dataset of adjacencies: simply scale them to [0, 1).
19

---

Published in Transactions on Machine Learning Research (06/2023)
In summary, the full procedure to produce [As+1)a; is as follows:
3 a A + B. (AoAk + AxAo) + 75.))
3 U jas ① (11 1)
force diagonal elements to 0 normalize entries per slice to be in [-1, 1]
U 3, U 3 max U J.::1
[A ReLU(U,,ar T/)
By normalizing in this way, we guarantee the intermediate matrix Un has entries in [-1,1] (before the ReLU). This plays two important roles. The first one has to do with training stability and to appreciate this point consider what could happen if no normalization is used. Suppose the entries of UjA,. are orders of magnitude larger than entries of U This can cause the model to push 7 >> Tf", leading to training instabilities and/or lack of convergence. The second point relates to interpretability of т. Indeed, the proposed normalization allows us to interpret the learned values 7 R on a similar scale. All the thresholds must be in [0, because: (i) anything above 1 will make the output all O; and (ii) we constrain it to be non-negative. In fact we can now plot all values (from all layers) against one another, and using the same scale ((0, 1)) interpret if a particular r is promoting a lot of sparsity in the output (r close to 1) or not (r close to 0), by examining its magnitude.
Gaining expressiveness by generalizing the channel reduction. MIMO filters provide a significant boost in performance, and effectively scale well to large networks due to their reliance on a minimal set of fast. operations such as batched (sparse) matrix multiplication and addition, as well as their extreme parameter efficiency. This is in contrast. to recent highly parametrized models e-g. Graph Transformers (Dwivedi et al.j 2022), which require significantly more computational resources effectively limiting their application to very small networks. Still, more expressiveness may be desired, and one such solution discussed below fits elegantly into the MIMO formulation of GDNs and still has nice scaling properties.
Consider a layer L which takes an input tensor Tin E R XNxN and outputs tensor Tout E X A single channel in the output filters T., to construct it's own filtered tensor of identical size. It then performs mean/sum reduction over the first. dimension of this constructed tensor, producing a single slice in the output tensor [i] N Note that this reduction essentially acts the filtered 'edge features' that the output channel produced. To gain expressiveness, this reduction transformation on the edge features can be swapped out for a more expressive transformation, e.g. an MLP, which is shared across all edges in the single output channel, and optionally all output channels within the layer (and optionally all layers in the network) reminiscent of the parameter sharing used in GNNs to transform node features. This sharing is what provides advantageous scaling properties.
Incorporating node features. In this work, we focus the network inverse problem introduced in Section 2 In this problem setting, we observe noisy graph Ao E NxN and look to recover the true latent graph Ar. E NRN, the problem does not consider the use of node features, and thus GDNs do not directly consider them. But if they are available, GDNs can readily leverage them with the use of MIMO filters and link-prediction for details. Given two nodes i.j E V with respective node features v, By E IR",, define a nonlinear function he : m x R. + 1Rm" with set of learnable parameters O As is standard in the Geometric Deep Learning literature [Bonsteiin et al. 2017], we can feed he any elementwise symmetric function of v, Uj c.g. the elementwise squared difference (v; u u,), elementwise Gaussian kernel of the difference g(v; vy), etc. We organize the pairwise outputs of he into a tensor R E R. XNxN such that R. R.j.a = l v,). In the MIMO formulation, the first. GDN layer already has a set of input channels Tin (some slices of T.n. can be prior information, learned, or zeros); we can simply append this new tensor R as more slices in Tin and proceed normally.
techniques which are standard in the DL literature; see the ensuing discussion
We leave these extensions for future work.
20

---

Published in Transactions on Machine Learning Research (06/2023)
A.5 Training
Training of the GDN model will be performed using stochastic (mini-batch) gradient descent to minimize task-dependent loss function _( as in (2). The loss is defined either as (i) the edgewise squared/ /absolute error between the predicted graph and the true graph for regression tasks, or (ii) a hinge loss with paramcter
2 0, both averaged over training set T == {A() A A ) namely
(0) ((A0 ()) = Σ ($(A (1) + (J - -) (-2(A0 )); 1)+
A =0
hingee (A
A Ann. >
masel A AL,0(A0 (8) 9)) 4 $(A0
1
maa (A 4 +)) 4 (0) (
I( T
Σ (
9 4(A 0)); u (hinge, mse, mae).
.ET
The loss is optimized with ADAM using learning rate of 0.01, 81 0.85, and B2 0.99.
Link prediction with GDNs and unbiased estimates of generalization. In the edge- regression lask, GDNs only use their validation data to determine when training has converged. When performing fink-prediction, GDNs have an additional use for this data: to choose the cutoff threshold E I2,, determining which raw outputs (which are continuous) should be considered positive odge predictions, at the end of training.
We use the training set to learn the parameters (via gradient descent) and to tune t. During the validation step, when then use this train-set-tunedl- on the validation data, giving an estimate of generalization error. This is then used for early-stopping, determining the best model learned after training, etc. We do not use the validation data to tune t during training. Only after training has completed, do we tune t with validation data. We train a handful of models this way, and the model which produces the best validation score (in this case lowest error) is the tested with the calidation-tunede thus providing an unbiased estimate of generalization.
A.6 Notes on the Experimental Setup
Synthetic graphs. For the experiments presented in Table I the synthetic graphs of size 68 are drawn from random graph models with the following parameters
Random geometric graphs (RG): 2, r 0.56.
Erdós-Rényi (ER): p p= .56.
Barabási-Albert (BA): 15
When sampling graphs to construct the datasets, we reject any samples which are not connected or have edge density outside of given range. For RG and ER, that range is [0.5,0.6], while BA the range is [0.3, 0.4]. This is an attempt to make the RG and ER graphs similar to the brain SC graphs, which have an average edge density of 0.56, and all SCs are in edge density range [0.5, 0.6). Due to the sampling procedure of BA, it is not possible to produce graph in this sparsity range, we lowered the range sightly. We thus take SCs to be an SBM-like ensemble and avoid a repetitive experiment with randomly drawn SBM graphs.
Note that the edge density ranges define the performance of the most naive of predictors: all ones/ ones/zeros. zeros. In the RG BG/BA/SC, an all ones prodictor achieves an average error of 44% (averra edge density). In the BAs, naive all zeros predictor achieves 35% -averrage edge density). This is useful to keep in mind when interpreting the results in Table 1
21

---

Published in Transactions on Machine Learning Research (06/2023)
Pseudo-syntheties. The 'pseudo-synthetic datasets' are those in which we diffuse synthetic signals over SCs from the HCP-YA dataset. This is an ideal setting to test the GDN models: we have weighted graphs to perform edge-regression on (the others are unweighted), while having Ao's that are true to our modeling assumptions. Note that SCs have a strong community-like structure, corresponding dominantly to the left and right hemispheres as well as subnetworks which have a high degree of connection, e-g. the Occipital Lobe which has 0.96 edge density - almost fully connected while the full brain network has edge density of 0.56.
Generation of the training set from synthetic diffused signals. We carry out the following steps to experiments Section 5.1 The
generate the training set := {40 " for our
:=
diffusion graph filter (:::) h)) fixed across all observations in the training set. We let K= and sample the filter coefficients h E I2 uniformly from the unit sphere. To generate each observation A we let P p= 50 and first sample w "...... ( ~ (O, IN) (i.i.d. standard Normal random vectors). For the given filter h, these white signals are diffused over A (using the graph models described above) to yield z/6 =H(A (0) h)w,) (6) for each p I....,,. P. We finally form the sample covariance
r Σ
A =
(2)00 p
(2)))
2 p p=1
This process is repeated for each sample i = 1,..., T in the training set. To examine robustness to the realizations of h, we repeat this data generation process three times (resampling the filter coefficients).
Normalization. Unless otherwise stated we divide Ao by its maximum eigenvalue before processing. A is used in its raw form when it is an unweighted graph. The SCs define weighted graphs and we scale all edge weights to be between 0 and by dividing by 9.9 (an upper bound on the maximum edge weight in the HCP-YA YA dataset). Similar to (Shrivastava al. 2020 Pu al. 2021 we find GLAD to be very sensitive to the conditioning of the labels AL; see the Appendix
A.7 for full discussion the normalization used in
this case.
Error and MSE of the baselines in Table 1 The edge weights returned by the baselines can be very small/large in magnitude and perform poorly when used directly in the regression task. We thus also provide a learned scaling parameter, tuned during the hyperparameter search, which provides approximately an order of magnitude improvement in MSE in GLASSO and halved the MSE Spectral Templates and Network Deconvolution. In link-prediction, we also tune a hard thresholding parameter on top of each method to clean up noisy outputs from the baseline, only if it improved their performance (it does). For complete details on the baseline implementation; see the Appendix A.7 Note also that the MSE reported in Table 1 is per-edge squared error averaged over the test set, namely M7 [Teass] E
[ (6)) ( where M == NNN - 1)
is the number of edges in a fully connected graph (without solf loops) with N = 68 nodes.
Size generalization. Something lo note is that we do not lune the threshold (found during training on the small N 68 graphs) on the larger graphs. We go straight from training to testing on the larger domains. Tuning the threshold using validation set (of larger graphs) would represent an casier problem. The model at no point, or in any way, is introduced to the data in the larger size domains for any form of maining/tuning
We decide to use the covariance matrix in this experiment, as opposed to the sample covariance matrix, as our Ao's 's This is for the simple reason that it would be difficult to control the SNR with respect to generalization error and would be secondary to the main thrust of the experiment. When run with number of signals proportional to graph size N, we see quite a similar trend, but due to time constraints, these results are not presented herein, but is apt for follow up work.
HCP data. HCP-YA provides up to 4 resting state fMRI scanning sessions for each subject, each lasting 15 minutes. We use the fMRI data which has been processed by the minimal processing pipeline (Van Essen et al. 2013). For every subject, after pre-processing the time series data to be zero mean, we concatenate all available time-series logether and compute the sample covariance matrix (which is what as used are Ao in the brain data experiment 5.2).
Due to expected homogeneity in SCs across a healthy population, information about the SC in the test set could be leveraged from the SCs in the training set. In plainer terms, the average SC in the training set is
22

---

Published in Transactions on Machine Learning Research (06/2023)
Parietal
Frontal
Occipital
Temporal
Unspecified
Figure 6: The four lobes in the brain cortex: parictal, frontal, occipital and temporal.
an effective predictor of SCs in the test set. In our experiments, we take random split of the 1063 SCs into 913/50/100 aaaiing/validationnttoss sets, and report how much we improve upon this predictor.
Raw data available from https://www.hummaanctt www.humanconneettone..rg.. org.//udy/ haauudy/hpppp hho--young--aull /overview, and subject to the HCP Data Use Agreement See here for discussion on the risks of personal identifiability.
Co-location and social networks among high school students. The Thiers13 dataset (Génois & Barral 2018) followed students in a French high school in 2013 for days, recording their interactions based on physical proximity (co-location) using wearable sensors as well as investigating social networks between students via survey. The co-location study traced 327 students, while a subset of such students filled out the surveys (156). We thus only consider the 156 students who have participated in both studies. The co-location data is a sequence of time intervals, each interval lists which students another in such interval. To construct the co-location network, we let vertices represent students and assign the weight an edge between student i and student to represent the (possibly 0) number of limes over the 5 day period the two such students were in close proximity with one another. The social network is the Faccbook graph between students: vertex i and vertex j have an unweighted edge between them the social network student and student j have are friends in Facebook. Thus we now have a co-location network Ao,ttal and social network Ar,total- To construct dataset of graph pairs = [Ao (1) A., (0) 7000 we draw random subsets of N 120 vertices (students). For each vertex subset, we construct a single graph pair by restricting the vertex set in each network to the sampled vertices, and removing all edges which attach to any node not in the subset. If either of the resulting graph pairs are not connected, the sample is not included.
within physical proximity to
We report performance of each model with an additional tuned threshold on the output, only if it increases performance. ND does not change the output mcaningfully (and thus the performance is identical to simply thresholding). GLAD continually converged to local minimum with worse performance than simple thresholding. Due the scaling issucs and resulting large computational costs a subsct of this datasct is used for ND (train/wa baain//al/aee /aal/val//ett test of 500 100/100/200 and GLASSO/Spel Temp (train/val/ten of 20/20/85) Increasing dataset sizes further did not seem to improve performance for any model.
A.7 Baselines
In the broad context of network topology identification, recent DL-based latent graph inference approaches
such DGCNN (Wang et. al. 2019), DGM (aazi et. al. 2023 NRI (Kipf al. 2018). PGN (Veličković
23

---

Published in Transactions on Machine Learning Research (06/2023)
et al. 2020 have been shown effective in obtaining better task-driven representations of relational data for machine learning applications, or to learn interactions among coupled dynamical systems. However, because the proposed GDN layer does not operate over node features, none of these state-of-the-art methods are appropriate for tackling the novel network doconvolution problem we are dealing with here. Hence, for the numerical evaluation of GDNs we chose the most relevant baseline models that we outline Section P and further describe the sequel in no particular order.
We were fair possible in the comparison of GDNs to baseline models. All baselines were optimized to minimize generalization error, which is what is presented in Table ! Many baseline methods aim to predict sparse graphs on their own, yet many fail to bring edge values fully to zero. We thus provide a threshold, tuned for generalization error using a validation set, on top of each method only if it improved the performance of the method in link- mre-prrdiition. The edge weights returned by the baselines can be very small/larg /large in magnitude and perform poorly when used directly the the regression task. We thus also provide a scaling parameter, tuned during the hyperparameter search, which provides approximately an order of magnitude improvement in MSE for GLASSO and halved the MSE for Spectral Templates and Network Deconvolution.
Unrolled PDS (L2G). L2G [Pu et al. 2021 is an unrolling of the dual splitting algorithm (PDS) introduced by (Kalofolia]] 2016). which aims to recover sparse, connected graphs from observations of smooth signals w.r.t. the latent graph At. More formally, consider the matrix [ #p] E NxP whose columns Ep are the observations in X. The rows, denoted by a. E R collect all P measurements at vertex i. Define then the nodal Fuclidean-distance matrix E E NNN where Egy- | x 11 V. Using these notions, and taking Ao E the signal smoothness measure of X over A can be equivalently written
primal-
1 Σ - Tr LX) = 2 1 @ Aollı
(12)
where L = diag(A1) A is the graph Laplacian and O denotes the elementwise product (Kalofolias) 2016 We can use the link between signal smoothness and edge sparsity in 12 to pose convex inverse problem
| @ Aoll - al log( (A1) + 9 | AEA
À E arg max
(13)
2
where a, 8 are tumable regularization parameters and A denotes the convex set of allowable adjacency matrices (hollow diagonal, non-negative entries, symmetric).
We implemented L2G following the released code from the authors Pu et al. 2021 which uses depth of 20, a learning rate of 0.01, a learning rate decay of 0.95, and the ADAM optimizer. L2G uses intermediate outputs in the loss, and discount these intermediate losses using a loss-discounting factor = 0.9 which we also use in our implementation. We tried both with and without intermediate losses and did not find a significant difference. We do not use the additional Variational Autoencoder introduced in (Pu et al. 2021 Finally, for the link-prediction task, a threshold is chosen with validation set to map the output to binary decision over edges in the same manner as used in GDNs and GLAD see Appendix A.5
Graph Convolutional Network applied to degrees of Ao (GCN-D). GCNs are simple graph neural
network models introduced in (Kipf & Welling 2017). which take node features X E RN a graph A E RNN, and uses weight matrices Wo, W1 to produce latent node representation Z = AReLU(AXX, W1. The graph À defined by the 'renormalization trick' Á := D- 1 ÄD 2, where À = A+ + and D diag (1). In our experiments we take our graph to be the observed graph A Ao, and our node signals to be the nodal degrees x Aol d (C =1), produeing the encoding model:
z = (o (Aod.))W
(14)
To decode the latent representations into edge predictions we use a simple outer product docoder A, = ZZT In (Kipf & Welling 2016) a sigmoid is used on the output of the inner product to produce a probability, which is then fed into binary cross entr loss function. We found including this sigmoid did not produce
24

---

Published in Transactions on Machine Learning Research (06/2023)
good graph estimates stable training. So we remove it and directly use the outer product in our loss, which for link-prediction was a hinge loss and for regression it was a MSE loss.
We use as close to identical architectural and training setup to (Kipf & Welling 2016) as possible: hidden /latent channel dimensions of 32/16, learning rate of 0.01, ADAM optimizer, training for 200 epochs (we use SGD with large batches- 20% of training data per batch instead of full gradient descent as in (Kipf & Welling, 2016) to work with our GPU). We found that the standard ReLU did not allow rich gradient information to flow to the first layers weights, and so we replace it with LeakyReLU with slope = 0.2 which increased performance. Additionally, the degree node signals d Aol result in very large values on many of our datasets, and so we |d]la A hyperparameter search was performed on the learning rate, but we could not find one which improved performance over the 0.01 used in (Kipf & Welling 2016). We performed further hyperparameter searches the channel widths (hidden, latent) E {(8, (16, 8), (32, 16), (64,32)), and found all worked similarly well, and so used (32, 16) as in (Kipf & Welling/ 2016). We also tested whether one-hot encoded features (included with and without the degree feature) could aid performance, but they did not.
normalize via dnorm= = d
Graph Isomorphism Network applied to degrees of Ao (GIN-D). GINs are graph ncural network models introduced in (Xu et al. 2019), which take node features X E RNx, a graph A E RVV and updates node representations as:
MLP() (11 + (8)))) +
Σ uEN(())
z ))
(15)
In the first layer, the node representation are set to the input node features. Here the input node signals to be æ = Aol= = d (C=1). Note our graphs are weighted, and so in the sum we also use the edge weights occurring in the graph. Self-loops can occur, for example as variances on the diagonal in constructed covariance matrices, but are not included in the neighborhood of a node. Similar to the GCN-D baseline, to decode the latent representations into edge predictions we use an outer product decoder Å. = ZzT, where Z is the matrix of latent node representations 2D)) in the final layer. For link prediction we additionally feed this through a sigmoid and then a binary cross entropy loss, for edge regression we feed this straight to an MSE loss.
We use close to identical architectural and training setup to (Xu et al. 2019) as possible: D == 5 GIN layers (including the imput layer) are applied, all MLPs have 2 layers, hidden channel dimensions of 64, the (()) in all layers are not learned and set to 0- - denoted GIN-O in (Xu et al. 2019). batch normalization is applied on each hidden layer, a learning rate of 0.01, ADAM optimizer, decay learning rate by 0.5 every 50 epochs, and training was performed for 350 epochs and stopped carly after 40 epochs if no improvement in loss was found. Additionally, the degree node signals d = Aol result in very large values on many of our datasets, and so we normalize via dnorm= = Tdlla"
T
Hard Thresholding (Threshold). The hard thresholding model consists of a single parameter T, and generates graph predictions as follows
Å. = I =11])
where 1(-} is an indicator function, and > denotes entry-wise inequality. For the synthetic experiments carried out in Section 5.1 to learn the structure of signals generated via network diffusion, Ao is either a covariance matrix or a correlation matrix. We tried both choices in our experiments, and reported the one that performed best in Table 1
Graphical Lasso (GLASSO). GLASSO is an approach for Gaussian graphical model selection (Yuan & Lin 2007 Friedman et al. 2008). In the context of the first application domain in Section 3 we will henceforth assume a zero-mean graph signal æ ~ NO, Σz). The goal is to estimate (conditional independence) graph structure encoded in the entries of the precision matrix θ = Σz1 To this end, given an empirical covariance matrix Ao == Σn estimated from observed signal realizations, GLASSO regularizes the maximun-likelihood & norm, yielding the convex problem
estimator of θ,, with the
e E arg max { det det e tracc () -ai|@|l1 al|@]l1].
(16)
0
25

---

Published in Transactions on Machine Learning Research (06/2023)
We found that taking the entry-wise absolute value of the GLASSO estimator improved its performance, and so we include that in the model before passing it through hard-thresholding operator
A. - -1(161=-11
One has to tune the hyperparameters a and T for link-prediction (and a third, the scaling parameter described below, for edgo-weight regression).
We used the sklearn GLASSO implementation found in: https://scikkit-eeann o/tpa///eeeemmo///aaaaaa//////// Bonerated/skleane wwvariinccee.ooo..o.. covariaanceeraapiccallll.a graphical Laalll_ass__b__ html
It is important to note that we do not use the typical rross-validation procedure seen with GLASSO. Typically, GLASSO is used in unsupervised applications with only one graph being predicted from P observations. In our application, we are predicting many graphs, each with P observations. Thus the typical procedure of choosing a using the log-likelihood [the non-regularized part of the GLASSO objective in (16)) over splils of the observed signals, not splils of the training sel, results in worse performance (and a different. a for cach graph). This is not surprising: exposing the training procedure to labeled data allows optimize for generalization. We are judging the models on their ability to generalize lo unseen graphs, and thus the typical procedure would provide an unfair advantage to model. While tried both sample covariance and sample correlation matrices as Ao, we found that we needed the normalization that the sample correlation provides, along with an additional scaling by the maximum magnitude eigenvalue, in order lo achieve numerical stability. GLASSO can take a good amount of time to run, and so we limited the validaiton and test set sizes to an even 100/100 split. With only 2 to 3 hyperparameters to tune, we found this was sufficient (no significant differences between validation and test performance in all runs, and when testing on larger graph sizes, no difference in generalization performance).
Unrolled AM GLASSO Objective (GLAD). GLAD is unrolling of an alternating minimization (AM) iterative procedure to form the GLASSO estimator (16); refer to (Shrivastava et al. 2020) for full details. GLAD inherits sensitivity to the conditioning of the labels from the constraint that θ be positive semi-definite in (6). In the experiments presented in (Shrivastava et al. 2020), precision matrices are sampled directly to be used as labels, and their minimum eigenvalue is set to 1 by adding an appropriately scaled identity matrix (diagonal loading). This minimum eigenvalue correction was attempted on the raw adjacency matrices A but did not result in convergent training on any data domain in our experiments. The best rosults (reported in Table [ were found corresponding to the minimum cigenvalue corrected Laplacian matrix representation of the latent graph, i.e., L, == diag( Ar.1) Ar. + I. The Laplacian L,, referred to as O in (Shrivastava et al. 2020) is used only in the the loss function to make training converge. The reported metrics in Table 1 disrogard the diagonal clements and compute error and MSE based on the off-diagonal entries, thus ensuring an equivalent comparison with all other methods. Finally, for the fink-prediction task, a threshold 7 is chosen with a validation set to map the output to binary decision over edges in the same manner as used in GDNN; see also Appendix A.5
We use the GLAD configuration with learning rate Ir = 0.1 and 8 = 0.9, 82= 0.999, the Pmn has 4 hidden layers, and Aan has 2 hidden layers. An intermediate loss discount factor = 0.8 worked reasonably well in the domains tested. ER graphs caused GLAD to converge to a trivial all zeros solution when D = 30. We increased the depth to D == 50 and performed a search over the learning rate and minimum eigenvalue m. (previously set to 1), resulting in a GLAD model (D = 50, br = -05, m. = 100, all other identical) that decrcased crror (by 1%) in the link-prediction task relative to the naive all zero output. MSE got worse in this configuration relative to the previous. In Table [ we report respective metrics from their best configurations: the error with the deeper configuration and MSE with the shallower configuration.
presented in (Shrivastava et al. 2020 == 30 layers, ADAM optimizer
When training on the SCs in Table I we incorporate prior information in GLAD by setting the prior to | = (A1) - A + I, where A is the adjacency matrix of the edgewise mean graph over all SCs in the training split. In other words we use the minimum cigenvalue corrected Laplacian matrix (as above) representation of A. GDNs simply use A as discussed in Section 3 the prior is (() = (Ao + I)-1, but found inferior performance.
We also tried the default setup, where
26

---

Published in Transactions on Machine Learning Research (06/2023)
GLAD in general has an O(D run time complexity, where D is the depth of GLAD, in both the forward and the backward pass, due to the computation of the matrix square root and its corresponding gradient. These can be approximated using e.g. Denman-Beaver iterations or Newton-Schulz iterations, which reduce the time complexity to O(DTN²), where T is the number of iterations run. These iterative methods still lack guarantees of convergence in general. We used such approximations Lions without nolicing reductions performance to make the experiments run in reasonable amounts of time. Even with such approximations GLAD has significant memory and runtime limitations. This results from the large D, moderate sized T (required for reasonable/staie stable approximations), as well the use of shared intermediate MLPs Pan and Ann the former of which is called N times (one for each entry in O) per layer. This forces batch sizes smaller than what would be desired and limits applicability on larger graph sizes. Note that in contrast the dominant term in the memory & runtime complexity of GDN layers are the (lypically sparse) matrix multiplications, making training and deployment on larger graph sizes realizable; see Figure 3
Network Deconvolution (ND). The ND approach is "a general method for inferring direct effects from an observed correlation matrix containing both direct and indirect effects" (Feizi et al. 2013). ND follows three steps: linear scaling to ensure all eigenvalues A of Ao fall in the interval A E [-11, 11, digen-decomposition of X
the scaled Ao ddiag(A) V-1, and deconvolution by applying f(A,) = 1-,, to all cigenvalues. We then construct our prediction as Å == V diagg (-)) In FFizzi et al., 2013), it is recommended a Pearson correlation matrix be constructed, which we followed. We applied an extra hard thresholding on the output, tuned for best generalization error, further increase performance. For each result shown 500 graphs were used in the hyperparameter search and 500 were used for testing.
Spectral Templates (SpecTemp). The SpecTemp method consists of two-step process whereby one: (i) first leverages the model (1) to estimate the graph eigenvectors V from those of Ao (the eigenvectors of A and Ao coincide); and (ii) combine V with a priori information about G (here sparsity) and feasibility constraints on A to obtain the optimal eigenvalues À of A = Vdiag(A)VT.
The second step entails solving the convex optimization problem
(c)
argmin ||||, s. to ||A Vdiag(A) <E 11 (A.))
(17)
We first perform a binary search on E R2 over the interval [0, 2] to find 6min, which is the smallest value With emin in hand, we now run an iteratively (t henceforth denotes
which allows a feasible solution to 17 iterations) re-weighted fn-oorm minimization problem with the aim of further pushing small edge weights to 0
711" (thus refining the graph estimate). Defining the weight matrix W |
E <R. R.N.N where 7, E R, are appropriately chosen positive constants and Að = A (emiin), we solve a sequence t = I,..., T of weighted L n-morrm minimization problems
A == argmin [ @ Al) S. to | Vdiag(A) 1 Emins Al 1, A {A.)
(18)
If any of these problems is infeasible, then A is returned where I {0, .,, obtained solution.
10. T) is the last successfully
Finally, a threshold T is chosen with a validation set to map the output to binary decision over edges, namely
A. = =
We solve the aforementioned optimization problems using MOSEK solvers in CVXPY. Solving such convex optimization problems can very computationally expensive/ /sow,, and because there are only two hyperpa- rameters, 100 graphs were used in the validation set, and 100 in the test set.
Least Squares plus Non-convex Optimization (LSOpt). LSOpt consists of first estimating the polynomial coefficients in ( via least-squares (LS) regression, and then using the found coefficients in the optimization of (3) to recover the graph A, Due to the collincarity between higher order matrix powers, ridge regression can be used in the estimation of the polynomial coefficients. If the true order K of the polynomial is not known ahead of time, one can allow h E R and add an additional form regularization
27

---

Published in Transactions on Machine Learning Research (06/2023)
Table 3: Runtime and memory usage analysis. GDN, 12G, and GLAD report time (s) for a forward and backward pass on single sample (batch size 32), as well as peak GPU memory consumption, using a T4 GPU on g4dn.xlarge AWS instance. SpecTemp reports inference time (s) for a single sample, and peak memory consumption, using 96 Core x86 Processor on ce.metal AWS instance. Results are taken after 40 epochs/samples to allow the running process to stabilize.
MODELS
TIME (8) PER SAMPLE
PEAK MEMORY USAGE (GB)
GDN-1
8.78-12122
1.28
GDN-8
3. 1.-3-32
12.84
L2G
2.28-3=as3
1.11
GLAD
6.78-3=ss3
11.73
SPECTEMP
11.3E-0ase-
3.84
for the sake of model-order selection. With h in hand, we optimize (3 non-convex problem due to the higher-order powers of A, using ADAM with a learning rate of 0.01 and a validation set of size 50 to tune X.
We start with the most favorable setup, using the true h in the optimization of (3), skipping the LS estimation step. Even in such a favorable setup the results were not in general competitive with GDN predictions, and so instead of further degrading performance by estimating h, we report these optimistic estimates of generalization crror.
A.8 Additional Numerical Experiments
Runtime and memory usage. In Section 5.1 we discuss the runtime complexity of GDNs and T.2G (quadratic) well as the runtime complexity of the cubic baselines (GLAD, GLASSO, ND, SpecTTemp). We further discuss in Section A.7 the approximations made to the matrix square root in GLAD to reduce its runtime complexity to be quadratic. Here we explore the actual run time (s) per sample and peak memory usage (for the full batch) for a handful of models taken from the runs on the = 68 RG graphs in Table D Note that GDNs, L2G, and GLAD are DL models trained with stochastic gradient descent, whereas SpecTemp is the solution to a sequence of optimization problems. We thus report the time to run the forward and backward pass over a batch of size 32 (we normalize by batch size to get a per sample time) for the former and time for the full inference over single sample in the latter. GDN-1 refers to GDN model with D = 8 and C²= 1, iee no MIMO channels, while GDN-8 refers to GDN model D = 8 and C = 8 the same model used for the experiments in Table [ GLAD, L2G, and SpecTTemp also follow the same setup used in Table [ see Section A.7 for details. The fundamental operations in GDNs are batched matrix multiplications and additions both of which are highly parallelizable and, importantly, have optimized implementations available to us. This is particularly true for the case of sparse batches of matrices (a setting where we can approach a linear runtime complexity). While we were able to reduce GLAD's runtime complexity to be quadratic, it still requires multiple internal MLPs and deeper unrolled network for good performance. SpecTemp relies on the full eigendecompositioon and solving multiple internal convex optimization ion problems each time it performs inference for sample- posing issues for real time applications or large networks.
The results in Table 3 show GDNs have the fastest runtimes, with GDN-1's being ~ 3z (8x) faster and using ~ 10x (10z) less peak GPU memory than GDN-8 (GLAD). GDN-8 about twice as fast as GLADD, with similar poak GPU memory usage. SpecTemp orders of magnitude slower and uses significantly more memory per sample- owing the cubic operations contains.
GDN hyperparameter search. To explore the effect of the dominant hyperparameters in GDNs, we perform a simple grid scarch over Depth (D) f3. 8, 20} and Channels per Layer (C) E (1, 8, 15) on RG graphs used in Table [ (train/val] val/test of 91 913/500/5000 found around the value of .01 there was no major effect on performance. As can be seen in Table | for each fixed depth, there monotonic improvement in performance (dccrcase in error) with an increasing number of channels. This occurs until our GPU runs out of memory (OOM). For each fixed channel width, there seems to be improvement in performance when depth is increased from D to D = 8, then proceeding
We did not include learning rale in the search as we
28

---

Published in Transactions on Machine Learning Research (06/2023)
Table Hyperparameter Search over Depth (D) and Channels per Layer (C) of a GDN network. Reported are the mean and standard error of the test crror (%) on the link prodiction task. OOM indicates "Out of Memory:
1 c
7.0E- 2=1sa
4.88-247e4
3.5E- 3.5
OOM
8
5.
22064 2.9E-213m4
2.2E-2
OOM
20
4.0
OOM
OOM
degradation (or lack of improvement) in performance when D 20. We thus choose C' 8 and D = 8 to get the best performance, while staying within the memory constraints of our hardware. These results would suggest further improvements in performance are readily attainable with more compute resources.
The full list of hyperparameters are follows: Learning rate, depth D, channel widths C, normalization of input Ao, normalization of intermediate outputs A(k]. In searching over these hyperparameters, values for the learning rate were & {le 6, le 4, le 2, le 1) D {3, 8, 20), C' {1,3, 8, 15), normalization for Ao, A[k] E (maximum eigenvalue, frobenius norm, maximum absolute value norm, 99th percentile norm). Optimal values are learning rate = 1e 2, D = 8, 8, and Ao, [kk normalization's 's of maximum eigenvalue, and maximum absolute value norm, respectively.
The procedure used find such hyperparameter were as follows. Initially fix depth D 10 with channel width C = 1. This the 10's of layers with good performance, while not approaching the super deep networks which can have gradient vanishing problems. A hyperparameter search over appropriate normalization Ao {maximum eigenvalue, frobenium norm, maximum absolute value norm, 99th percentile norm} was performed first. Maximum eigenvalue normalization was found to perform best and resulted in the most stable training. Next, a search over appropriate normalization for intermediate outputs A[&] was performed over the same normalizations. We found both maximum eigenvalue and maximum absolute value norm perform similarly, but maximum eigenvalue normalization has higher computational cost and maximum absolute value norm projected the edges onto the interval [-1, 1] which was convenient for edge prediction applications. Next a search was performed over learning rate E { 1e 6, le - 4, le 2, le 1). D E {3, 8, 20} and C' e {1,3, 15). This search grid is displayed in Table 4 in A.8 with the learning rate dimension omitted for case of interpretation. All scarches were porformed over RG graphs from Table 1
inspired by unrollings in other domains which tend to use in
A.9 Source Code and Configuration
Code is available on github at https://github. cwwp:/awwwwa////p((SI We rely on PyTorch Paszke et al. 2019) (BSD license) heavily and use Conda (Anaconda] 2020) (BSD license) to make our system as hardware- independent as possible. GPUs were used to train the larger GDNs, which are available to use for free on Google Colab. Development was done on Mac and Limux (Ubuntu) systems, and not tested on Windows machines. The processed brain data was too large to include with the code, but is available upon request.
29


Reviews:
[{"review_id": "w0akeS7Srn", "reviewer": "Reviewer_VMWN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "1. This paper is an empirical work, but the introduction lacks motivation/intuition. More concretely, the first two paragraphs list the algorithms or conclusions of the existing works, e.g., the graph structure can be further optimized for downstream tasks. However, the use cases of such algorithms are not well discussed. What are the real applications and corresponding impacts that we need to optimize the downstream graph structure? And how to interpret the modification to the graph structure? Adding such explanations can help the audience know better about the importance of this task. Currently, this information is only vaguely discussed in Introduction.\n\n2. Also, in the Introduction, in the second paragraph, the authors mention they propose GDN, but what is the relation between GDN and related works? This is not explained in the Introduction.\n\n3. In Sec 2, I’m wondering if the authors can explain what is the relation between E and A_O and A_L? Now the formulation is not straightforward without a detailed explanation/intuition.\n- According to the context, I guess that E is equivalent to A_O, and they are known as inputs; A_L is the latent variable that we want to learn. Is this correct?\n- Also, can authors also help explain what are the real scenarios corresponding to A_O and A_L? First, the draft says A_O contains “spurious, indirect connections”, but such a description is not very straightforward for readers. Then, the authors list some use cases in the last paragraph in Sec 2, but they are not self-contained and reader-friendly to audiences not familiar with the domain-specific field. I would recommend authors explain this into more detail, and a simple case can be something like what are the A_O and A_L for each application.\n- The authors have some examples in Sec 3, I’m wondering what the relationship is between the examples in Sec 2 and 3.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper studies the graph structure learning/recovery problem and it can have the potential to solve many real applications (more explanations are required). Given the core graph decomposition equation in Eq (1), the authors propose a Graph Deconvolution Network (GDN) algorithm. GDN utilizes an iterative optimization manner for learning and unrolling for inference.\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. The writing can be improved. There is some key information missing in the draft, and some contents seem to be redundant (or more explanations on the difference are required). Please check the comments above.\n\n2. The authors can polish the Introduction section, especially on the motivation and technical novelty of GDN. Current writing cannot reveal the importance of the problem and contributions of this work.\n\n3. The equation references are missing. It should be “Eq/Equation (x)” instead of “(x)”.\n", "broader_impact_concerns": "none", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "bzYXlUvtqe", "reviewer": "Reviewer_AruC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths.**\n\n- [St1] The proposed method is derived by a principle idea of approximating the optimization process ([So1]).\n- [St2] The proposed method is generalized to graphs whose size is different from those in the training data set ([So3, N2]).\n- [St3] In the experiments, the proposed method is compared with many baselines, and the experimental conditions of the baselines are taken into account so that the evaluation is fair. ([So4])\n\n**Weaknesses**\n\n- [W1] Writing of related work has room for improvement ([C1, C2]) \n- [W2] No comparison with existing algorithms for solving link prediction tasks ([So5])\n- [W3] There is room for improvement in the numerical experiments to justify the linearization in the proposed method ([So2]).\n\n\n**Soundness**\n\n[So1] The derivation of GDN is theoretically sound and reasonable. That is, the architecture of GDN comes from estimating the gradient of the loss function by a single layer of GNN.\n\n[So2] The ablation study (p.11) wrote that $K\\geq 2$ resulted in unstable training (p.11. Section 5.1). I would like to clarify how the authors judged this.\n\n[So3] This paper claimed the advantages of the proposed method as follows:\n1. GDN can incorporate prior knowledge by giving an initial guess of the true graph as an initial value.\n2. By using GNNs, GDN can make inference on graphs whose size is different from those used during training.\n\nThe first point was verified by the ablation study (P10), although I think this is a general property of L2O methods ([N3]).\nThe second point was verified in Figure 3. Indeed, it has been claimed that the increase in error is modest with respect to graph size. Intuitively, I agree with this claim. However, it would be desirable to provide objective metrics to judge whether this performance degradation is really modest, such as comparison with baselines.\n\n[So4] Numerical experiments have compared GDN with many baselines. The experimental settings of the baselines are carefully designed to make the comparison fair. In addition, since the baseline methods are described in detail, their reproducibility is high.\n\n[So5] The proposed method has not been compared with existing methods for estimating graph structure from node features, such as those mentioned on p.9, because it was difficult to compare with them (Section 5 Baselines). However, I wonder if we can do so in the task of estimating the true graph structure from node features. For the proposed method, we should construct an empirical variance-covariance matrix or the precision matrix and solve the problem of estimating the true graph structure from it, similar to the graph diffusion model.\n\n**Clarity**\n\n[C1] Most of the paper was clear. There was no great difficulty in understanding it. However, I think there is room for improvement in Related work in section 1 and Remark in section 3.\n\n[C2] Although various concepts and existing methods were explained in the related work section, I felt that they were not well organized,\n- The relationship between existing and proposed methods was only implicitly stated. For example, in the last part of the related work, some studies on latent graph learning and algorithm unrolling were mentioned. But we could not see why they were there in the first reading.\n- Deep generative models have been mentioned. However, I think they are less relevant to the proposed method because they are not used in the supervised setting nor in graph recovery problems.\n\n**Novelty and Significance**\n\n[N1] In this paper, the authors considered the problem of estimating the true graph structure from an observed graph. This problem setting includes the problems considered in Pu et al. (2021) and Shrivastava et al. (2019). To the best of my knowledge, there are no known papers that apply GNNs to these problems. In this sense, this paper is novel.\n\n[N2] The proposed method allows generalization to test data of different sizes, which was not possible with previous techniques. Therefore, I think this point is significant, although I think it would be desirable if we could verify it in a more objective way.\n\n[N3] The fact that GDN can introduce prior knowledge as an initial guess has been mentioned as an advantage. However, I think this property is common to other L2O methods, although they are not applicable to the problem considered in this paper. Therefore, the novelty of the proposed method is limited from this point of view.\n\n[N4] The technique of replacing the gradient of the objective function with a linearized model to reduce the computational complexity was a new technique in the context of L2O. On the other hand, as pointed out in [So2], it would be desirable to improve the validation of its effectiveness.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper considered supervised graph learning problems where the observed graph is represented by a polynomial of the true graph structure.\nFor this problem, the paper proposed Graph Deconvolution Network (GDN). This model is obtained by linearizing a single step of optimization by proximal gradient descent of an L1-regularized loss function.\nSeveral architectural improvements (prior optimization, multi-input multi-output filter, and decoupled layer paramters) were also proposed to improve the accuracy.\nThe proposed method was applied to the tasks of link prediction and edge weight regression on synthetic data and two types of real data (neuroimaging data and social network data), and the usefulness of the proposed method was verified by conducting an ablation study.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "[R1] I would like to clarify what the authors meant by \"the novel supervised setting for topology identification\" in the Summary of Main Contributions was the content of the Remark in Section 3. In other words, this paper considered the problem of estimating the true graph from the observed graph, which included the problem of estimating the graphical model from an empirical variance-covariance matrix (or precision matrix).\n\n[R2] That paper claimed that \"the GDN learns the distribution of the graph\" (e.g., in the abstract). I would like to clarify that this claim meant that the GDN learns $P(y|x)$, where $x$ is the corrupted graph and $y$ is the unknown true graph structure.\n\n[R3] Similar terms have been used in different contexts, which may confuse readers who are not familiar with this area. Specifically, I would like to clarify whether \"graph structure learning\" (abstract), \"network topology inference\" (P2), \"latent graph learning\" (P2), \"topology identification\" (P2), and \"graph recovery\" (abstract) all referred to the same task or different concepts.\n\n[R4] P10, Section 5.1: I would like to clarify how the initial value of the block diagonal initialization was determined.\n\n[R5] Section 5.1: I would like to clarify how to generate $A^{(i)}\\_{O}$ (i.e. $\\hat{\\Sigma}\\_{x}^{(i)}$) from $A^{(i)}\\_{L}$. Is the following correct?\n\n$w^{(i)}_1, \\ldots, w^{(i)}_P \\sim N(0, I)$\n\n$x^{(i)}_p = H(A^{(i)}L)w^{(i)}_p$\n\n$\\hat{\\Sigma}\\_{x}^{(i)} = \\frac{1}{P} \\sum\\_{p=1}^P x^{(i)}\\_p x^{(i)\\top}\\_p$\n\n[R6] Some citations were not formatted correctly. For example, Bronstein et al. (2017) and Ortega et al. (2018) in Section 1 should be replaced with (Bronstein et al. ) and (Ortega et al., 2018), respectively (i.e., use citep).\n\n[R7] I think we can simplify the notation $\\bm{H}(\\bm{A}; \\bm{h})$ to $\\bm{H}(\\bm{A})$ because the second argument $\\bm{h}$ has not changed in this paper. We can interpret it as an assignment of $A$ to $\\bm{H}(x) = \\sum_{k=0}^K h_k x^k$.\n\n[R8] The E notation in Table 1 and Table 2 was different.\n\n[R9] References: Some references did not give the source (e.g., Pu et al. (2021), Kipf & Welling (2016)).", "broader_impact_concerns": "This paper does not have major concerns about broader impacts.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "qjVP3zpxt3", "reviewer": "Reviewer_XTUC", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* The investigation of the network inverse/deconvolution problem is interesting, which is claimed to have broad applications across various domains, ranging from inferring protein contact structure to epidemiology.\n* The discussion on the problem background is comprehensive, which is very helpful for justifying the applicability of GDNs. The two discussed properties of GDNs, i.e., symmetry preserving and node permutation equivariance, are closely connected with the family of GNNs.\n* The paper is generally well-written and almost clear everywhere. I enjoyed reading through this manuscript.\n\nWeaknesses:\n* The outline of the manuscript is a bit weird given the related work is embedded in the Introduction Section, but this does not really hurt the understandability.\n\n* One concern from mine is the naming of network deconvolution, which is not very clearly defined in my opinion. In other literature such as [1], the process of deconvolution on graphs is explicitly denoted as recovering the original node features/signals given smoothed node representations obtained by a GNN. This deconvolution on the node features/signals is different from the notation of network deconvolution in this manuscript, which is indeed a network reverse problem and focuses on the graph structures. I suggest the authors define the task of network deconvolution more explicitly, with possible examples from the literature. The discussion on this difference would help audiences not be confused about the task dealt with in this paper.\n\n* Since the architecture of GDNs does not involve the node features, could the authors provide some insights on how the node features on graphs might affect the performance of GDNs? For example, if there exist two graphs with different node features but the same network structure $\\boldsymbol A_{O}$\nreveal different latent graph $\\boldsymbol A_{L}$ , could GDNs distinguish the two graphs during either the training or inference stage?\n\n* For the empirical evaluation, the discussion on runtime and memory usage is favored. I am curious about why L2G[2] is not included in either complexity analysis or numerical comparison since L2G and Glad are the most related attempts to GDNs.\n\n[1] Deconvolutional Networks on Graph Data, NeurIPS 2021\n\n[2] Learning to learn graph topologies, NeurIPS 2021\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper systematically studies the problem of how to infer the graph structure from (partially) observed graphs as an optimization problem, which they denote as a network inverse/deconvolution problem. To solve this problem, the authors propose a neural network-based solution, GDN, from the principle of algorithm unrolling, and use truncated proximal projection. GDNs are differentiable with respect to the neural network parameters and the input graphs and are proven to have a node permutation equivariance property. Empirical evaluations on both synthetic and two real-world datasets, i.e., neuroimaging and social network datasets, demonstrate the effectiveness of GDNs on both link prediction and edge-weight regression tasks.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "My concerns are mainly from the mighty misleading name deconvolution, the discussion on how node features affect the performance of GDNs, and the runtime and memory usage comparison with L2G. Please kindly refer to Section Strengths and Weaknesses for detailed suggestions.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the supervised graph structure learning/recovery problem. More specifically, the paper considers the problem of learning an adjacency matrix from a noisy input adjacency matrix A_O.  Then, the paper proposed Graph Deconvolution Network (GDN). This model is obtained by linearizing a single step of optimization by proximal gradient descent of an L1-regularized loss function. 

Overall, the idea is solid and well supported by the numerical experiments; it is good to accept to TMLR. 
Please prepare for a camera-ready version by including the author's name and updating the character colors.
Decision: Accept (poster)