A Model-Agnostic Graph Neural Network for Integrating Local and Global
  Information
Graph Neural Networks (GNNs) have achieved promising performance in a variety
of graph-focused tasks. Despite their success, however, existing GNNs suffer
from two significant limitations: a lack of interpretability in results due to
their black-box nature, and an inability to learn representations of varying
orders. To tackle these issues, we propose a novel
\textbf{M}odel-\textbf{a}gnostic \textbf{G}raph Neural \textbf{Net}work
(MaGNet) framework, which is able to effectively integrate information of
various orders, extract knowledge from high-order neighbors, and provide
meaningful and interpretable results by identifying influential compact graph
structures. In particular, MaGNet consists of two components: an estimation
model for the latent representation of complex relationships under graph
topology, and an interpretation model that identifies influential nodes, edges,
and node features. Theoretically, we establish the generalization error bound
for MaGNet via empirical Rademacher complexity, and demonstrate its power to
represent layer-wise neighborhood mixing. We conduct comprehensive numerical
studies using simulated data to demonstrate the superior performance of MaGNet
in comparison to several state-of-the-art alternatives. Furthermore, we apply
MaGNet to a real-world case study aimed at extracting task-critical information
from brain activity data, thereby highlighting its effectiveness in advancing
scientific research.
Wenzhuo Zhou, Annie Qu, Keiland W. Cooper, Norbert Fortin, Babak Shahbaba
# A Model-Agnostic Graph Neural Network for Integrating Local and Global Information

###### Abstract

Graph Neural Networks (GNNs) have achieved promising performance in a variety of graph-focused tasks. Despite their success, existing GNNs suffer from two significant limitations: a lack of interpretability in results due to their black-box nature, and an inability to learn representations of varying orders. To tackle these issues, we propose a novel **M**odel-**a**gnostic **G**raph Neural **N**etwork (MaGNet) framework, which is able to sequentially integrate information of various orders, extract knowledge from high-order neighbors, and provide meaningful and interpretable results by identifying influential compact graph structures. In particular, MaGNet consists of two components: an estimation model for the latent representation of complex relationships under graph topology, and an interpretation model that identifies influential nodes, edges, and important node features. Theoretically, we establish the generalization error bound for MaGNet via empirical Rademacher complexity, and showcase its power to represent layer-wise neighborhood mixing. We conduct comprehensive numerical studies using simulated data to demonstrate the superior performance of MaGNet in comparison to several state-of-the-art alternatives. Furthermore, we apply MaGNet to a real-world case study aimed at extracting task-critical information from brain activity data, thereby highlighting its effectiveness in advancing scientific research.

**Keywords:** Graph representation; Empirical Rademacher complexity; Information aggregation

## 1 Introduction

Graph-structured data is ubiquitous throughout the natural and social sciences, from brain networks to social relationships. In the most general view, a graph is simply a collection of nodes representing entities such as people, genes, and brain regions, along with a set of edges representing interactions between pairs of nodes. By representing such interconnected entities as
graphs, it is possible to leverage their geometric topology to study functional relationships in systems with network based frameworks. For example, when studying brain function, graphs can be used to model activity relationships among neurons or brain regions. For such applications, it is essential to build relational inductive biases into machine learning methods in order to develop systems that can learn, reason, and generalize. To achieve this goal, in recent years there has been a surge in research on graph representation learning, including techniques for deep graph embedding powered by deep learning architecture, e.g., neural message passing approaches inspired by belief propagation. These advances in graph representation learning have led to new developments in various domains, including the analysis of information processing in the brain, chemical synthesis, recommender systems, and modeling social networks.

Among the methods for deep graph representation, the family of graph neural networks (GNNs) has achieved remarkable success in real-world graph-focused tasks (Defferrard et al., 2016; Zhang and Chen, 2018). In general, GNNs iteratively aggregate and combine node representations within a graph, so-called message passing, to generate a set of learned hidden representation features (Gilmer et al., 2017). The main neural architectures of GNNs include graph convolutional networks (GCNs; Kipf and Welling (2016)), graph attention networks (GATs; Velickovic et al. (2017)), and graph isomorphism networks (GINs; Xu et al. (2018a)), and many others variants (Hamilton, 2020).

While GNNs are capable of capturing subgraph information through local operations, e.g., graph convolution operations, they can be prone to over-smoothing the learned representations when applying multiple rounds of local operators. This can lead to models that treat all nodes uniformly, resulting in node representations that converge towards an indistinguishable vector (Li et al., 2018). Another issue with over-smoothing is that it limits the ability to capture high-order information, which can be only aggregated by executing a sufficient number of local operations, such as when nodes have a large receptive field. Several studies have suggested that this over-smoothing phenomenon is a key factor in the performance degradation of deep GNNs (Bodnar et al., 2022). Therefore, most recent works (Zhang et al., 2022; Keriven, 2022) advocate the use of various shallow GNNs (e.g., up to three layers) to avoid the over-smoothing issue.
Unfortunately, as previously discussed, simply shallowing GNNs often falls short in capturing high-order information due to insufficient rounds of message passing.

Another key factor for learning sufficient representation power in GNNs involves directly incorporating and sequentially combining information from neighbors at different node neighbor orders, including both low-order and high-order levels. Statistically, by incorporating both low-order (immediate neighbors) and high-order (i.e., neighbors beyond the immediate vicinity) information, GNNs can learn a richer representation under graph topology. For example, in the brain, low-order information may capture relationships within a functional region and high-order information relationships across functional regions (Bassett and Sporns, 2017). Importantly, the information flows in a sequential manner starting from within-region functionality to cross-region functionality (Vazquez-Rodriguez et al., 2019).

Models that adhere to the principle of effectively combining different-order information are known as multi-scale GNNs, as they enable the exploration and integration of information at multiple levels of granularity within a graph (Xu et al., 2018; Liao et al., 2019; Sun et al., 2019; Oono and Suzuki, 2020; Liu et al., 2022). Multi-scale GNNs provide enhanced flexibility and adaptability in facilitating information flow across various neighborhood ranges. This is achieved by directing the outputs of intermediate layers to contribute to the final representation. Nonetheless, existing methods struggle to effectively integrate representations of different orders in a sequential manner due to their memoryless property: the most recent learned representation tends to overwrite the previously acquired representation. In addition, this approach still suffers from over-smoothing issues and fails to capture high-order latent representations in general.

To address these issues, we propose a novel **M**odel-**a**gnostic **G**raph neural **Net**work (MaGNet) framework with two major components: the _estimation model_ and the _interpretation model_. The estimation model captures the complex relationship between the feature information and a target outcome under graph topology, allowing for powerful latent representation. The interpretation model, on the other hand, identifies a compact subgraph structure, i.e., influential nodes and edges and a small subset of node features that have a crucial role in the learned estimation model. This provides a meaningful explanation of the underlying mechanisms of the relationships,
allowing for a better understanding of the factors governing the behavior of complex dynamic systems. For example, for brain activity data, the interpretation model may be used to identify the most informative patterns of regional activity and of functional relationships among spatially distributed regions, as well as the most informative time periods.

The main advantages of the proposed framework and our contributions are summarized as follows. First, the proposed neural architecture of the model alleviates the over-smoothing issue and thus can effectively extract knowledge from high-order neighbors. In addition, the designed actor-critic neural architecture has the ability to integrate various order information via resolving the memoryless issue. It sequentially combines the representation learned from actor graph neural networks, while a critic neural network plays a role in evaluating the quality of the hidden representation learned by the actor network.

Second, we propose a tractable and concise interpretation framework for the black-box estimation model. We formulate our interpretation model as an optimization task that maximizes the information gained between a prediction of the estimation model and the distribution of possible subgraph structures. The approach is model-agnostic and can explain the estimation model on any graph-focused task for graphs without assumptions regarding the underlying true model or data-generating mechanisms.

Third, we study the ability to integrate various-order information as well as the statistical complexity of the proposed model via a minimax measurement on empirical Rademacher complexity. Unlike existing analyses solely working on standard message passing neural networks, our results can be applied to a mixture of neural architectures between message passing neural networks and feedforward neural networks, where the message passing neural network is a special case in our theoretical framework. Furthermore, we provide a rigorous generalization error bound for the MaGNet estimation model that consists of sequential deep-learning component models.

The rest of the paper is organized as follows. Section 2 lays out the basic graph notation and model settings. Section 3 formally defines the estimation module. Section 4 illustrates how to construct the interpretation model for an estimated model. A comprehensive theoretical study of the proposed framework is provided in Section 5. Sections 6 and 7 demonstrate the empirical
performance of our methods. We conclude this paper with a discussion of possible future research directions. All technical proofs are provided in Appendix.

## 2 Preliminaries

### Graph Structure

In this section, we present some preliminaries and notations used throughout the paper. Let \(G=(V,E)\) represent the graph, where \(V\) represents the vertex set consisting of nodes \(\{v_{1},v_{2},...,v_{N}\}\), and \(E\in V\times V\) denotes the edge set with \((i,j)\)-th element \(e_{ij}\). The number of total nodes in the graph is denoted by \(N\). A graph can be described by a symmetric (typically sparse) adjacency matrix \(A\in\{0,1\}^{N\times N}\) derived from \(V\) and \(E\). In this setting, \(a_{ij}=0\) indicates that the edge \(e_{ij}\) is missing, whereas \(a_{ij}=1\) indicates that the corresponding edge exists. There is a \(T\)-dimensional set of features, \(X_{i}\), associated with each node, \(v_{i}\) so that the entire feature set is denoted as \(X\in\mathbb{R}^{N\times T}\). Suppose we have observed \(n\) graph instances, each consisting of a fixed graph structure but with different node features. Let \(G_{i}\) denote the \(i\)th instance of a graph, where \(i\in{1,2,\ldots,n}\). While our approach can be used for predictive models in general, here we focus on classification problems, where the objective is to assign a binary label \(s\in\{-1,1\}\) to each graph instance.

### Neural Message Passing

The basic graph neural network (GNN) model can be motivated in a variety of ways. The same fundamental GNN model has been derived as a generalization of convolutions to non-Euclidean data (Bodnar et al., 2022), and as a differentiable variant of belief propagation (Dabkowski and Gal, 2017), as well as by analogy to classic graph isomorphism tests (Graham et al., 2019). Regardless of the motivation, the defining feature of a GNN is that it uses a form of neural message passing in which vector messages are exchanged between nodes and updated using neural networks (Abu-El-Haija et al., 2019). During each message passing iteration in a GNN, a hidden embedding corresponding to each node \(v\in\mathcal{V}\), denoted as \(H_{v}^{(k)}\), is updated according to information aggregated from \(v\)'s graph neighborhood \(\mathcal{N}(v)\). This message passing update can be
expressed as follows:

\[H_{v}^{(k+1)} =f_{\text{update}}^{(k)}\left(H_{v}^{(k)},f_{\text{agg}}^{(k)}\left( \left\{H_{u}^{(k)},\forall u\in\mathcal{N}(v)\right\}\right)\right)\] \[=f_{\text{update}}^{(k)}\left(H_{v}^{(k)},M_{\mathcal{N}(v)}^{(k) }\right),\]

where \(f_{\text{update}}\) and \(f_{\text{agg}}\) are the update and aggregate functions, which are arbitrary differentiable functions (i.e., neural networks). The term \(M_{\mathcal{N}(v)}\) is the "message" that is aggregated from \(v\)'s graph neighborhood \(\mathcal{N}(v)\). We use superscripts to distinguish the embeddings and functions at different iterations of message passing. At each iteration \(k\) of the GNN, the aggregate function takes as input the set of embeddings of the nodes in \(v\)'s graph neighborhood \(\mathcal{N}(v)\) and generates a message \(M_{\mathcal{N}(v)}^{(k)}\) based on this aggregated neighborhood information. The update function then combines the message \(M_{\mathcal{N}(v)}^{(k)}\) with the previous embedding \(H_{v}^{(k-1)}\) of node \(v\) to generate the updated embedding \(H_{v}^{(k)}\).

### 2.3 Graph Convolutional Network

Let \(D\) be the degree matrix corresponding to the adjacency matrix \(\tilde{A}=A+I\) with \(D_{ii}=\sum_{j=1}^{N}\tilde{A}_{ij}\). The hidden graph representation of nodes with two graph convolutional layers (Kipf and Welling, 2016) can be formulated in a matrix form:

\[H=\widetilde{\mathcal{L}}\text{ReLU}(\widetilde{\mathcal{L}}XW^{(0)})W^{(1)}, \tag{1}\]

where \(H\in\mathbb{R}^{N\times T^{(1)}}\) is the final embedding matrix of nodes and \(T_{hidden}\) is the dimension of the node hidden representation. The graph Laplacian is defined as \(\widetilde{\mathcal{L}}=D^{-\frac{1}{2}}\tilde{A}D^{-\frac{1}{2}}\). In addition, the weight matrix \(W^{(0)}\in\mathbb{R}^{t_{latent}\times t^{(0)}}\) is the input-to-hidden weight matrix for a hidden layer with \(t^{(0)}\) feature maps, and \(W^{(1)}\in\mathbb{R}^{t^{(0)}\times t^{(1)}}\) is the hidden-to-output weight matrix. Here we consider the two-layer case that aims to simplify the notation; the above definition can be easily extended to \(l\) graph convolutional layers with \(l>2\).
Estimation Model

In this section, we introduce a novel graph neural network, which aims to represent feature information and explain an outcome of interest by effectively analyzing the relationships between them. To achieve this goal, it integrates both low-order and high-order neighbor node information in order to form a powerful latent representation. Here, the low-order information refers to information aggregated from the local neighbors of a node, while the high-order information means the messages aggregated beyond just the immediate neighbors, capturing the global graph information.

To characterize the feasibility of capturing various-order information rigorously, we first introduce a generalized second-order \(\Delta\)-operator which was defined in Abu-El-Haija et al. (2019), i.e., \(K\)-order counterparts for \(K\geq 3\) as follows:

**Definition 3.1**.: _Given a graph neural network, \(\Delta(K)\) representer represents \(K\)-order node neighbor information for \(K\in\mathbb{N}\), e.g., there exists a real-valued vector \(\nu=(\nu_{1},\nu_{2},...,\nu_{K})\) and an injective (one-to-one) mapping function \(g(\cdot)\), such that the output embedding of this graph neural network can be represented as_

\[g\left(\sum_{k=0}^{K}\nu_{k}\cdot\mathcal{L}^{k}X\right):=\Delta(K),\]

_for any type of graph Laplacian \(\mathcal{L}\) operation and input node feature matrix \(X\)._

Learning such an operator enables GNNs to represent feature differences among \(K\)-order node neighbor's information. When a candidate GNN model learns the \(\Delta(K)\)-operator, it effectively captures the \(K\)-order order neighborhood information in the hidden representation.

In GCNs, the graph representation is obtained through interactions of neighboring nodes during multiple rounds of learned message passing. Ideally, one could consider a deep architecture via stacking \(K\) GCN layers in order to learn a \(\Delta(K)\) operator. Despite this, most of the GCN models employ shallow architectures, typically utilizing only two or three-order information (Zhang et al., 2020). The reason behind this limitation is two-fold. First, when repeatedly applying Laplacian smoothing, GCNs may mix node features from different clusters, rendering
them indistinguishable. This phenomenon is known as the _over-smoothing_ issue (Li et al., 2021). Second, most GCNs are built upon a feedforward mechanism and suffer from the _memoryless_ problem. After each layer operation, the representation learned from the current layer overwrites the representation produced from the previous layers, meaning that there's no explicit memory mechanism. In other words, the _over-smoothing_ issue creats difficulties in capturing high-order information, while the _memoryless_ issue leads to a loss of lower-order information. Theoretically, under Definition 3.1, the GCN models with over-smoothing or memoryless issues cannot learn the \(\Delta(K)\)-operator in the case where \(K\) is large.

### Actor-Critic Graph Neural Network

In this section, we propose an actor-critic graph neural network that is designed to aggregate different levels of node-neighbor information to compute a powerful graph embedding. In this dual neural network structure, the actor graph neural network aims to capture the hidden representation for each order of node-neighbor information, while the critic neural network plays a role in evaluating the quality of the hidden representation learned by the actor network. With the hidden representation and evaluated quality scores, we perform a fusion operation to integrate the representation using calculated quality scores as weights. This framework resolves the over-smoothing and memoryless issues, which is guaranteed to learn a \(\Delta(K)\)-operator.

In contrast to GCNs, we adopt the simple weighted sum aggregator and abandon the use of feature transformation and nonlinear activation. As a result, the graph convolution operation in our actor graph neural network is defined as:

\[H^{(k)}=(\mathcal{L})^{k}XW, \tag{2}\]

where the graph Laplacian \(\mathcal{L}=D^{-\frac{1}{2}}AD^{-\frac{1}{2}}\) and the weight matrix \(W\) is fixed to be an identity matrix in every message passing round. We argue that the nonlinear feature transformation is not critical and that the majority of the benefit arises from the local averaging of neighboring features. This is because the vectorized temporal signal is not like the multi-dimensional image data that require the use of many nonlinear layers to capture sufficient information. By removing the contracting mapping non-linear activation, our graph convolution layer achieves slower
convergence in certain values for embedding vectors, thus alleviating the over-smoothing issue. Additionally, the abandoned feature transformation operation greatly facilitates computation. It is also worth noting that in (2), we aggregate only the connected neighbors and do not integrate the target node itself (i.e., self-connection). That is, the graph Laplacian is based on the adjacency matrix \(A\) instead of its augmented counterpart \(\widetilde{A}\). This is different from most existing graph convolutions (Kipf and Welling, 2016) that typically aggregate extended neighbors and need to handle the self-connection specially. The fusion operation, to be discussed later, essentially captures the same effect as self-connections. Further, based on the node representation \(H^{(k)}\), we can apply a graph pooling operation to summarize the graph embedding from \(H^{(k)}\). In the following, we use the notation \(H^{(k)}\) for graph embedding.

As previously stated, our primary objective is to aggregate mixed-order information. To this end, we propose a fusion operation to completely preserve the information from different-order neighbors. Specifically, we can regard the graph embedding \(H\) from \(k\)-th round of message

Figure 1: An illustrating example of 3-layer neural architecture of the actor-critic graph neural network.
passing as the output of the \(l\)-th order information summarization, denoted as \(H^{(k)}\). From the perspective of meta algorithms, the embedding \(H^{(l)}\) can be regarded as the graph embedding learned from the \(k\)-th actor graph neural network. Then the ultimate graph embedding can be weighted combined from the various order graph embeddings, i.e.,

\[\widetilde{H}=\sum_{k=1}^{K}\alpha^{(k)}H^{(k)}, \tag{3}\]

where \(\alpha^{(k)}\) is the fusion weight corresponding to the quality (or importance) of the \(k\)-th order knowledge for \(k=1,...,K\). This fusion operation can be understood as the ensemble of multiple single actor networks, i.e., the actor networks for generating graph embedding \(H^{(k)}\). Thus, the fusion operation naturally combines the unique characteristics of different single learners such as the different order information.

To determine the fusion weights \(\alpha^{(k)}\), we introduce a critic network, e.g., feedforward neural network, which plays a role in evaluating the quality of \(H^{(k)}\) in a bias induction way. That is

\[\alpha^{(k)}=\frac{1}{2}\log\left(\frac{1-\epsilon^{(k)}}{\epsilon^{(k)}} \right),\]

where the function \(\text{logit}(x)=\log(x/(1-x))\) and the error rate \(\epsilon^{(k)}\) is defined as

\[\epsilon^{(k)}=\sum_{i=1}^{n}\beta_{i}^{(k)}\mathds{1}\left\{s_{i}\neq\underset {\{m=-1,1\}}{\arg\max}\text{softmax}(H_{i}^{(k)})\right\}/\sum_{i=1}^{n}\beta _{i}^{(k)}. \tag{4}\]

Here \(H_{i}^{(k)}\) denotes the graph embedding for the \(i\)-th graph sample, \(\arg\max_{\{m=-1,1\}}(x)\) is the operator taking the maximum element in the two-dimensional vector \(x\), and the coefficient \(\beta_{i}^{(k)}\) denotes the weight of the \(i\)-th graph sample. Intuitively, the \(\epsilon^{(k)}\) can be understood as the weighted classification error rate for the \(k\)-th single actor network, i.e. \(\text{softmax}(H^{(k)})\). For \(i=1,...n\), we adjust the graph sample weights \(\beta_{i}^{(k)}\) sequentially from the \(k\)-th to the \((k+1)\)-th step by following the updating rule:

\[\beta_{i}^{(k+1)}\propto\beta_{i}^{(k)}\exp\left(\alpha^{(k)}\mathds{1}\left\{ s_{i}\neq\underset{\{m=-1,1\}}{\arg\max}\text{softmax}(H_{i}^{(k)})\right\} \right).\]

This update rule intentionally pays more attention to the misclassified graph samples with potentially insufficient representation power, and increases their weights when training the next
single actor network.

Looking at the ultimate graph embedding \(\widetilde{H}\) in (3), we can see that it adaptively combines information from the first to the \(K\)-th order node neighbors information, as opposed to simply mixing different order embedding with equal weights. Furthermore, the sequential updates of single actor networks maintain the same learning pattern as standard GCNs, wherein the message passing for the \((k+1)\)-th hop directly succeeds the message passing for the \(k\)-th hop. Due to this, our actor-critic graph neural network preserves most of the desirable properties in standard GCNs, such as invariance to graph isomorphism (Xu et al., 2018) and relational representation (Wu et al., 2020).

Ultimately, we derive the classification model and prediction rule given the \(i\)-th graph embedding \(\widetilde{H}_{i}\),

\[p(\cdot|\widetilde{H}_{i})=\text{softmax}(\widetilde{H}_{i}). \tag{5}\]

This equation returns the classification logits. To streamline the notation throughout the paper, we use \(p_{\widehat{\theta}}(\cdot)\) to represent a trained actor-critic graph neural network model. To have a better understanding of the proposed actor-critic neural network, we illustrate its neural architecture in Figure 1.

Finally, we note that the model training in (5) is a well-studied convex optimization problem. It can be performed using efficient second-order methods or stochastic gradient descent (SGD) (Bottou, 2010). As long as the graph connectivity pattern remains sufficiently sparse, SGD can naturally scale to handle very large graph sizes. Furthermore, we ensure the neural network architecture's consistency by training the layers sequentially, following the order of node neighbors' message passing. This sequential training approach enables us to utilize the trained parameters from the previous training to initialize the current model training, which effectively reduces computational costs.
Interpretation Model

Although the estimation model boasts representation power to capture the complex relationship between the outcome of interests and features, understanding the rationale behind its predictions can be still quite challenging. In this section, we present a novel interpretation framework designed to uncover the reasoning behind the "black-box" estimation model.

To bridge the gap between estimation and interpretation, we first observe that our estimation model extracts temporal feature information from various-order node neighbors as well as graph topology to output the hidden graph representation for predictions. This suggests that the prediction made by the estimation model, i.e., \(\widehat{s}=\arg\max_{\{m=-1,1\}}p_{\widehat{\theta}}(\cdot)\) in (5), is determined by the adjacency matrix \(A\) and the node feature information \(X\). Formally, to comprehend the model mechanism and provide explanations, the problem is transformed into identifying important/influential subgraphs, denoted as \(G_{sub}\subseteq G\) with a corresponding adjacency matrix \(A_{sub}\), and a small subset of the full-dimension of the node feature \(X_{sub}\). We first focus on the identification of influential subgraphs by assuming \(X_{sub}\) has been obtained, and then discuss how to perform node feature selection simultaneously with subgraph identifications.

We adapt the principle of information gain, which was first introduced in the context of decision trees (Larose and Larose, 2014), into our framework. In particular, we formulate an optimization framework for influential subgraph identification. Our goal is to maximize the information gain with respect to subgraph candidates \(G_{\text{sub}}\) in order to determine the most influential one:

\[\operatorname*{argmax}_{G_{sub}}\operatorname{IG}\left(p_{\widehat{\theta}}, G_{sub}\right)=\eta(p_{\widehat{\theta}})-\eta\left(p_{\widehat{\theta}}\mid G _{sub},X_{sub}\right), \tag{6}\]

where \(\eta(\cdot)\) and \(\eta(\cdot|\cdot)\) denote the entropy and conditional entropy respectively.

Essentially, information gain quantifies the change in prediction probability between the full model \(p_{\widehat{\theta}}(\cdot)\) and the one constrained to the subgraph \(G_{sub}\) and the subset node feature \(X_{sub}\). For example, if removing edge \(e_{ij}\), i.e., the \((i,j)\)-th element in the adjacency matrix \(A\), from the full graph \(G\) significantly decreases the prediction probability, then this edge is influential and should be included in the subgraph \(G_{sub}\). Conversely, if the edge \(e_{ij}\) is deemed redundant for prediction
by the learned estimation model, it should be excluded.

Examining the right-hand side of (6), we can easily observe that the entropy term \(\eta(p_{\widehat{\theta}})\) remains constant since the parameters \(\widehat{\theta}\) are fixed for an estimated model. Consequently, the objective of minimizing information gain in (6) is equivalent to maximizing the conditional entropy \(\eta\left(p_{\widehat{\theta}}\mid G_{sub},X_{sub}\right)\). Nevertheless, directly optimizing the above objective function is intractable, as there are \(2^{|V|}\) candidates for the subgraph \(G_{sub}\). To address this issue, we consider a relaxation by assuming that the subgraph is a Gilbert random graph (Gilbert, 1959). The selection of edges from the original input graph \(G\) are conditionally independent of each other and follow a probability distribution. In detail, the edge \(e_{ij}\) is a binary variable indicating whether the edge is selected, with \(e_{ij}=1\) if selected and 0 otherwise. Therefore, the graph \(G_{sub}\) is a random graph with probability

\[P(G_{sub})=\Pi_{i,j\in N}P(e_{ij}).\]

A straightforward instantiation of \(P(e_{ij})\) is the Bernoulli distribution \(e_{ij}\sim\text{Bern}(\mu_{ij})\), where \(\mu_{ij}\) is the first moment. In particular, we can rewrite the parametrized objective as:

\[\underset{G_{sub}}{\text{Minimize }}\eta\left(p_{\widehat{\theta}}\mid G_{ sub},X_{sub}\right)=\underset{G_{sub}(\mu)}{\text{Minimize }}\mathbb{E}_{G_{sub}(\mu)}[\eta\left(p_{\widehat{\theta}}\mid G_{ sub},X_{sub}\right)], \tag{7}\]

where \(G_{sub}(\mu)\) is the parametrized random subgraph. Due to the discrete nature of the subgraph \(G_{sub}(\mu)\), the objective function is non-smooth, making optimization challenging and unstable. To address this issue, we further propose a continuous approximation for the binary sampling process. Let \(\epsilon\) be a uniform random variable, i.e., \(\epsilon\sim\text{Unif}(0,1)\), and the real-valued parameters \(\psi_{ij}\in\Psi\), and a temperature parameter \(\omega\in\mathbb{R}^{+}\), then a sample of the binary edge \(e_{ij}\) can be approximated by a sigmoid mapping:

\[\widetilde{e}_{ij}=\text{sigmoid}\left(\frac{\log(\epsilon)-\log(1-\epsilon) +\psi_{ij}}{\omega}\right).\]

We denote \(\widetilde{G}_{sub}(\Psi)\) as the continuous relaxation counterpart of the subgraph, with the \((i,j)\)-th element of the adjacency matrix being \(\widetilde{e}_{ij}\). Interestingly, in our analysis, the temperature parameter \(\omega\) can describe the relationship between \(\widetilde{G}_{sub}(\Psi)\) and \(G_{sub}(\mu)\). We observe that as
\(\omega\to 0\), the approximated edge \(\widetilde{e}_{ij}\) converges to the edge \(e_{ij}\), with the probability mass function,

\[\lim_{\omega\to 0}P(\widetilde{e}_{ij}=1)=\frac{\exp(\psi_{ij})}{1+\exp(\psi_{ij})}.\]

Recall that the edge \(e_{ij}\) follows a Bernoulli distribution with mean \(\mu_{ij}\). If we reparameterize \(\psi_{ij}\) such that:

\[\psi_{ij}=\log\left(\frac{\mu_{ij}}{1-\mu_{ij}}\right),\]

we achieve asymptotical consistency of the approximated subgraph, i.e., \(\lim_{\omega\to 0}\widetilde{G}_{sub}(\Psi)=G_{sub}(\mu)\). This supports the feasibility of applying continuous relaxation to the binary distribution.

Unlike the objective function in (7) induced by the discrete original subgraph, the objective function becomes smooth under the edge continuous approximation and can be easily optimized using gradient-based methods. In other words, the gradient of the continuous edge approximation \(\widetilde{e}_{ij}\) with respect to the parameters \(\psi_{ij}\) is computable. More importantly, the sampling randomness towards the subgraph is absorbed into a uniform random variable \(\epsilon\) peel-off from the parameterized binary Bernoulli distribution, which greatly relaxes the complexity of the sampling processing.

In this manner, the objective function in (7) can be reformulated as

\[\underset{\Psi}{\text{Minimize}}\ \mathbb{E}_{\epsilon\sim\text{Unif}(0,1)}[ \eta\left(p_{\widehat{\theta}}\mid G_{sub}(\Psi),X_{sub}\right)].\]

Unfortunately, solving the conditional entropy is still computationally expensive. To avoid this issue, we follow Kipf et al. (2018) to minimize a cross-entropy as the objective function. The empirical objective then becomes

\[\underset{\Psi}{\text{Minimize}}\ \frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\sum_{m \in\{-1,1\}}p_{\widehat{\theta}}(s_{i}=m|X_{sub}(i))\log p_{\widehat{\theta} }(s_{i}=m|G_{sub}(\Psi),X_{sub}(i)),\]

where \(n_{0}\) is the sampling size and \(p_{\widehat{\theta}}(s_{i}=m|G_{sub}(\Psi),X_{sub}(i))\) denotes the classification logits conditional on the subgraph \(G_{sub}(\Psi)\) and the subset feature \(X_{sub}(i)\) of the \(i\)-th graph sample.

So far, in the above development, we implicitly assume that the subset feature \(X_{sub}\) is known, which is not the case in practice. In the context where the subset feature \(X_{sub}\) is not given, the main challenges are: 1) identifying the subset is unknown; and 2) the fact that integrating
this feature selection into the developed subgraph identification optimization framework is not trivial. Motivated by the great success of self-supervised techniques in large neural language models (Devlin et al., 2018), we propose to use a "masking" approach to convert the feature subsetting problem into an optimization problem that can be naturally combined with subgraph identification. Specifically, we define a binary vector \(\mathcal{B}\in\{0,1\}^{T}\) which holds the same dimension as the raw node feature. For each node \(v_{i}\) and its raw node feature \(X_{i}\), \(i=1,...,N\), we multiply the raw feature with the binary vector \(\mathcal{B}\), i.e., \(X_{i}\odot\mathcal{B}\), where \(\odot\) is the Hadamard product. Intuitively, the vector \(\mathcal{B}\) converts the value in some dimension of the node feature to 0. This aligns with the rationale that if a particular feature is not important, the corresponding weights in the neural network weight matrix take values close to 0. In terms of the principle of information gain, this type of masking does not significantly decrease the probability of the prediction or alter the information gain.

The binary vector \(\mathcal{B}\) is non-smooth; we also consider a continuous relaxation for the vector by leveraging a sigmoid mapping, so that the feature selection procedure becomes a smooth optimization problem, i.e.,

\[X\odot\text{sigmoid}(\widetilde{\mathcal{B}}),\]

where \(\widetilde{\mathcal{B}}\in\mathbb{R}^{T}\) is a real-valued vector and the \(\text{sigmoid}(\widetilde{\mathcal{B}})\) is applied to each row of \(X\). Next, we remove the low values in \(\widetilde{\mathcal{B}}\) through thresholding to arrive at the feature subsetting.

Now, the subgraph identification and the feature selection can be naturally integrated into a single minimization problem:

\[\min_{\Psi,\widetilde{\mathcal{B}}}\,\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\sum_{m \in\{-1,1\}}p_{\widetilde{\theta}}(s_{i}=m|X(i)\odot\text{sigmoid}(\widetilde {\mathcal{B}}))\log p_{\widetilde{\theta}}(s_{i}=m|G_{sub}(\Psi),X(i)\odot \text{sigmoid}(\widetilde{\mathcal{B}})),\]

which forms a unified optimization framework. This unified optimization framework allows for the simultaneous identification of influential subgraphs and important node features, leading to a more interpretable and efficient model. The resulting optimization problem can be solved using gradient-based techniques.
Theory

In this section, we present the main theoretical results of the estimation model. First, we demonstrate that in contrast to standard GCNs, the proposed estimation model is capable of representing the feature differences among various-order neighbors. Second, we study the capacity of the actor graph neural network in terms of empirical Rademacher complexity. The derived bound is minimax optimal through careful analysis of both the lower and upper bounds. Furthermore, we provide a probabilistic upper bound on the generalization error of the actor-critic graph neural network which is calibrated in the fusion algorithm.

**Theorem 5.1**.: _The MaGNet actor-critic graph neural network is capable of learning a \(\Delta(K)\)-representer, i.e., it is able to sufficiently capture \(K\)-order node neighbor information._

Theorem 5.1 demonstrates that our estimation model can learn various-order information. This ensures the capability of the proposed estimation model on high-order message passing, where nodes receive latent representations from their 1-order neighbors as well as further \(K\)-order neighbors at the information aggregation step. In contrast, the existing GCNs are not capable of representing this class of operations, even when stacked over multiple layers. To establish the bounds on generalization errors and simplify the proof, we make the following technical assumptions.

**Assumption 5.1**.: _The \(L_{2}\)-norm of the feature vector in the input space is bounded, namely, for some constant \(\widetilde{c}>0\), we assume \(\left\|X_{i}\right\|_{2}\leq\widetilde{c}\) for all \(i=1,...,N\)._

**Assumption 5.2**.: _The maximum number of elements in the graph Laplacian matrix is bounded above by, i.e., \(\max_{i\in[N]}\max_{j\in[N]}|\mathcal{L}_{ij}|\leq c_{\mathcal{L}}\)._

**Assumption 5.3**.: _The Frobenius norm of any weight matrix in the estimation model is bounded. Namely,_

\[\|W^{(l)}_{\text{MLP}}\|_{F}\leq c_{2},\|W^{(l_{0})}_{\text{MLP}} \|_{F}\leq c_{1},\|W\|_{F}\leq c_{0},\]

_with some constant \(c_{0},c_{1},c_{2}>0\) and the Frobenius norm \(\|\cdot\|_{F}\), where \(\|W^{l}_{\text{MLP}}\|_{F}\) is the weight matrix in \(l\)-th layer of MLP for any \(l<l_{0}\)._
**Assumption 5.4**.: _The element in the graph Laplacian \(\mathcal{L}\) is bounded by \(c_{2}>0\), namely, \(|\mathcal{L}_{ij}|\leq c_{2}\) for any \(i,j\)._

**Assumption 5.5**.: _The number of neighbors of each node is equal to each other, namely, for some common constant \(q\in\mathbb{N}^{+}\), assume \(q:=N(v_{i})\) for all node \(v_{i}\in V\), where \(N(\cdot)\) indicating the neighborhood nodes._

The above assumptions are common in the (graph) neural network literature. Assumption 5.1, 5.2 and 5.3 impose norm constraints on the parameters, graph laplacian matrix, and input feature, making the model class fall into a compact metric space (Liao et al., 2020). Assumption 5.4 is a standard assumption to control the intensity of the graph Laplacian in GNN literature (Hamilton, 2020; Lv, 2021). Assumption 5.5 requires us to focus on homogeneous graphs. We first present our result on bounding the Rademacher complexity of the model class \(\mathcal{F}_{c_{0},c_{1},c_{2}}\), which is the estimation model part before and up to the step that produces \(H^{(K)}\). For \(i_{0}\)-th graph sample, formally, we define our estimation model class \(\mathcal{F}_{c_{0},c_{1},c_{2}}\) in the setting of \(K=3\) and \(l_{0}=2\) without loss of generality,

\[\mathcal{F}_{c_{0},c_{1},c_{2}}:=\bigg{\{}f(X(i_{0}))=\sigma \Bigg{(}\sum_{q=1}^{d_{1}}w_{\mathrm{MLP}}^{(2)}q\sigma\Bigg{(}\sum_{t=1}^{k }w_{\mathrm{MLP}}^{(1)}t_{q}\frac{1}{N}\sum_{m=1}^{N}\sum_{i=1}^{N}\mathcal{L} _{mi}\sum_{v=1}^{N}\mathcal{L}_{iv}\] \[\qquad\qquad\qquad\qquad\times\sum_{j\in N(v)}\mathcal{L}_{vj} \left<X(i_{0})_{j},\mathbf{w}_{t}\right>\Bigg{)}\Bigg{)},\quad i_{0}\in[n],\] \[\|W_{\mathrm{MLP}}^{(1)}\|_{F}\leq c_{2},\|W_{\mathrm{MLP}}^{(2) }\|\leq c_{1},\|W\|_{F}\leq c_{0}\bigg{\}},\]

where \(\sigma(\cdot)\) is some identical mapping function, \(w_{\mathrm{MLP}}^{(1)}\) and \(w_{\mathrm{MLP}}^{(2)}\) is the first and second layer of the MLP for critic network, and \(\mathbf{w}_{t}\) is the \(t\)-th row of the weight matrix \(W\). Note that we use this particular setting as an example of the model class for simplifying the expression. The following theoretical results hold for the general case of \(K\) and \(l_{0}\).

**Definition 5.1**.: _Given the input node feature matrix \(\{X(i)\}_{i=1}^{n}\) and the \(K\)-layers actor-critic graph neural network class \(\mathcal{F}_{c_{0},c_{1},c_{2}}\), the empirical Rademacher complexity of \(\mathcal{F}\) is defined as_

\[\widehat{\mathcal{R}}(\mathcal{F}_{c_{0},c_{1},c_{2}}):=\mathbb{E}_{\epsilon} \left[\frac{1}{n}\sup_{f\in\mathcal{F}}\bigg{|}\sum_{j=1}^{n}\epsilon_{j}f \left(X(j)\right)\bigg{|}X(1),X(2),\ldots,X(n)\right],\]
_where \(\{\epsilon_{i}\}_{i=1}^{n}\) is an i.i.d. family of Rademacher variables, independent of \(\{X(i)\}_{i=1}^{n}\)._

**Theorem 5.2**.: _Under Assumptions 5.1-5.5, the empirical Rademacher complexity is bounded by_

\[\widehat{\mathcal{R}}(\mathcal{F}_{c_{0},c_{1},c_{2}})\leq \frac{(L_{0})^{l_{0}}c_{0}c_{1}c_{2}(\min_{i,j\in[N]}\mathcal{L}_ {i,j})^{K-1}\widetilde{c}q^{K}}{\sqrt{n}}|\lambda_{\min}(\mathcal{L})|; \text{Upper Bound}\] \[\widehat{\mathcal{R}}(\mathcal{F}_{c_{0},c_{1},c_{2}})\geq \frac{2^{K}(L_{0})^{l_{0}}c_{0}c_{1}c_{2}c_{\mathcal{L}}^{K-1} \widetilde{c}q^{K+0.5}}{\sqrt{n}}|\lambda_{\max}(\mathcal{L})|; \text{Lower Bound}\]

_where \(L_{0}\) is the Lipschitz constant for non-linear activation function in the critic network, \(\lambda_{\min}(\mathcal{L})\) and \(\lambda_{\max}(\mathcal{L})\) are the finite minimum and maximum absolute eigenvalue of graph Laplacian \(\mathcal{L}\) and \(c_{0},c_{\mathcal{L}},q,\widetilde{c}\) are defined in the empirical Rademacher complexity upper bound._

Theorem 5.2 demonstrates that our derived upper bound is tight up to some constants when comparing it to the lower bound. Theorem 5.2 indicates that the upper bound of \(\widehat{\mathcal{R}}(\mathcal{F}_{c_{0},c_{1},c_{2}})\) depends on the number of graph instances, the degree distribution of the graph, and the graph convolution filter. Interestingly, the above bound is independent of the maximum number of nodes \(N\), for traditional regular graphs.

Applying our results in empirical Rademacher complexity \(\widehat{\mathcal{R}}(\mathcal{F}_{c_{0},c_{1},c_{2}})\) to generalization analysis, we now state the fundamental result of the generalization bound of the estimation model. We denote \(conv(\mathcal{F}_{c_{0},c_{1},c_{2}})\) as the closed convex hull of \(\mathcal{F}_{c_{0},c_{1},c_{2}}\). That is, \(conv(\mathcal{F}_{c_{0},c_{1},c_{2}})\) consists of all functions that are pointwise limits of convex combinations of functions from \(\mathcal{F}\):

\[\text{conv}(\mathcal{F}_{c_{0},c_{1},c_{2}}):=\Big{\{}f:\forall x,f (x)=\lim_{K\rightarrow\infty}f_{K}(x),f_{K}=\sum_{k=1}^{K}w_{k}f_{k},\] \[\sum_{k=1}^{K}w_{k}=1,f_{k}\in\mathcal{F}_{c_{0},c_{1},c_{2}},K \geq 1\Big{\}}.\]

Obviously, we can observe that the combination in (3) belongs to \(conv(\mathcal{F}_{c_{0},c_{1},c_{2}})\). Next, we present the probabilistic generalization error for the estimation model:

**Theorem 5.3**.: _Under Assumptions 5.1-5.5, given the \(\widehat{s}\) is the predicted label from actor-critic graph neural network class with \(K\) layers and true label \(s\), then the probabilistic upper bounds of
the generalization error_

\[P\left(\widehat{s}s_{0}\leq 0\right)\leq\prod_{k=1}^{K}\Bigg{\{} \underbrace{2\sqrt{\epsilon^{(k)}\left(1-\epsilon^{(k)}\right)}+ \left(\frac{\log\log_{2}\left(2\big{(}\log\prod_{k=1}^{K}\sqrt{\frac{1- \epsilon^{(k)}}{\epsilon^{(k)}}}\lor 1\big{)}\right)}{n}\right)^{1/2}}_{\text{fusion estimation bias}}+\underbrace{\sqrt{\frac{1}{2n}\log\frac{2}{\delta}}}_{\text{ intrinsic uncertainty}}\] \[+\underbrace{\frac{2^{K+3}(L_{0})^{l_{0}}c_{0}c_{1}c_{2}c_{ \mathcal{L}}^{K-1}\widetilde{c}q^{K+0.5}}{\sqrt{n}}|\lambda_{\max}(\mathcal{L} )|\Bigg{(}\log\prod_{k=1}^{K}\sqrt{\frac{1-\epsilon^{(k)}}{\epsilon^{(k)}}} \lor 1\Bigg{)}}_{\text{local complexity}}\Bigg{\}},\]

_with probability at least \(1-\delta\) for \(\delta\in[0,1)\)._

Theorem 5.3 demonstrates that the generalization error of the proposed estimation model is bounded, in terms of the error rate at the \(k\)-th iteration defined in (4), i.e., \(\epsilon^{(k)}\geq 1/2\) for any \(k=1,...,K\). In comparison to the generalization bound on vanilla GNNs (Scarselli et al., 2018; Garg et al., 2020), our bound is independent of the number of hidden units and the maximum number of nodes \(N\) in any input graph. For a regular graph with \(q=\mathcal{O}(1)\)(Bollobas, 1998), we conclude that \(\lambda_{\max}(\mathcal{L})=1\), which yields a generalization error bound of order \(\mathcal{O}(1/\sqrt{n})\) that is fully independent of the number of nodes \(N\).

## 6 Simulation Studies

In this section, we present a comprehensive evaluation of the MaGNet using synthetic datasets. Specifically, we investigate its accuracy via binary classification tasks. Furthermore, we evaluate the performance of MaGNet's interpretation model by conducting thorough experiments for various purposes (e.g., node-wise, edge-wise, and feature-wise).

To generate graphs, we allow the number of nodes \(V\) in the graph to vary, ranging from 30, 50, to 75, with the graph sample sizes \(n\) being 100 or 250. Each node has a temporal feature of dimension \(p\) as 20 or 50. We distinguish two categories of nodes, specifically, important nodes and non-important nodes, and we generate their features by applying two separate processes, resulting in two different settings.

**Setting 1**: For the important nodes, features are generated following a multivariate Gaussian distribution, \(\text{MVN}(0,0.1\cdot I_{p\times p})\), where \(I_{p\times p}\) is an identity matrix. On the other hand, for the non-important nodes, the instance features are sampled from a uniform distribution, \(\text{Unif}(0,0.5)\)
**Setting 2**: For the important nodes, features are generated following a Gaussian process with the mean function as \(\mu(x_{t})=0.1(t-25)^{2}\) and its kernel covariance function as \(k\left(x_{t},x_{t+h}\right)=\exp\left(-\left|x_{t}-x_{t+h}\right|^{2}/(2\sigma^ {2})\right)\) to have a dependency between the temporal features, where \(\sigma=0.5\) and \(h=3\). On the other hand, for the non-important nodes, the instance features are sampled from a Gaussian process with the mean function as \(\mu^{\prime}(x_{t})=0.1t\) with the same kernel covariance function \(k\left(x_{t},x_{t+h}\right)\) as in the important nodes generation.

The distinct generation of features for each node class aims to create a distribution gap influencing the classification target outcome. It is important to note that the target outcome or classification rule is based solely on the features of important nodes and remains independent of those of non-important nodes. This generation process allows a good classifier to be able to separate important nodes from non-important ones. In all settings, we consider a binary outcome of interest, i.e., \(s=\{-1,1\}\), with the classification rule \(X_{V_{0}}/|V_{0}|+N(0,0.1)>0\), where \(N(0,0.25)\) introduces noise into the classification and \(X_{V_{0}}\) are the features of the important nodes and \(|\cdot|\) is the cardinality operator.

### Evaluation on Estimation Model

In this section, we evaluate the classification accuracy of the MaGNet estimation model, comparing it against several benchmark approaches, including DeepNet (LeCun et al., 2015), penalized-logistic regression (Hastie et al., 2009), and RandomForest (Breiman, 2001). Furthermore, we compare our method with Graph Neural Networks (GNNs) models such as the vanilla graph convolutional network (V-GCN) (Kipf and Welling, 2016), graph attention network (GAT) (Velickovic et al., 2017), and the Graph Isomorphism Network (GIN) (Xu et al., 2018a). To accommodate the competing methods for non-GNN methods, we transform the unstructured graph data into structured data. Specifically, we treat each node as a distinct sample within the dataset. Under this setting, we have a total of \(n|V|\) samples for training sets. This pre-processing ensures a fair comparison between our model and the non-GNN competing methods. The results are provided in Table 1 and Table 2.

As illustrated in Tables 1 and 2, the MaGNet estimation model is the best classifier among all the methods. This superiority is consistent across varying sample sizes, node quantities, and
important node sizes, indicating MaGNet's robust performance in binary classification tasks. The performance gains of the MaGNet model mainly come from the ability to integrate both local and global information, which yields powerful representations for graph-structured data.

### Evaluation on Interpretation Model

In this section, we first introduce three types of interpretation tasks. Then we report the performance of the MaGNet interpretation model with comparisons to the competing methods.

#### 6.2.1 Interpretation Tasks

As mentioned above, we consider three types of model interpretation tasks: node-wise, edge-wise, and feature-wise interpretation tasks. We note that each of them is aligned with the functionalities of the proposed MaGNet interpretation model. In the following interpretation tasks, we focus

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline
**Sample Size** & **Important Nodes** & **Nodes** & **MaGNet** & **V-GCN** & **GAT** & **GIN** & **DeepNets** & **RF** & \(l_{1}\)**-Logit** \\ \hline \multirow{2}{*}{100} & \multirow{2}{*}{10} & 30 & **0.779** & 0.755 & 0.762 & 0.758 & 0.732 & 0.721 & 0.685 \\  & & 50 & **0.757** & 0.727 & 0.747 & 0.741 & 0.721 & 0.704 & 0.665 \\  & & 75 & **0.742** & 0.712 & 0.722 & 0.724 & 0.703 & 0.688 & 0.651 \\ \cline{2-10}  & \multirow{2}{*}{20} & 30 & **0.791** & 0.775 & 0.768 & 0.764 & 0.757 & 0.752 & 0.704 \\  & & 50 & **0.768** & 0.748 & 0.754 & 0.737 & 0.714 & 0.701 & 0.678 \\  & & 75 & **0.758** & 0.738 & 0.730 & 0.729 & 0.711 & 0.693 & 0.667 \\ \hline \multirow{2}{*}{250} & \multirow{2}{*}{10} & 30 & **0.788** & 0.758 & 0.766 & 0.769 & 0.739 & 0.722 & 0.694 \\  & & 50 & **0.767** & 0.733 & 0.741 & 0.727 & 0.710 & 0.707 & 0.661 \\  & & 75 & **0.759** & 0.744 & 0.745 & 0.732 & 0.713 & 0.700 & 0.668 \\ \cline{2-10}  & \multirow{2}{*}{20} & 30 & **0.797** & 0.767 & 0.754 & 0.748 & 0.728 & 0.710 & 0.684 \\  & & 50 & **0.771** & 0.741 & 0.756 & 0.733 & 0.752 & 0.747 & 0.697 \\  & & 75 & **0.771** & 0.750 & 0.745 & 0.727 & 0.733 & 0.718 & 0.687 \\ \hline \hline \end{tabular}
\end{table}
Table 1: The results of classification accuracy over 50 repeated experiments in Setting 1.

\begin{table}
\begin{tabular}{c c c c c c c c c c} \hline \hline
**Sample Size** & **Important Nodes** & **Nodes** & **MaGNet** & **V-GCN** & **GAT** & **GIN** & **DeepNets** & **RF** & \(l_{1}\)**-Logit** \\ \hline \multirow{2}{*}{100} & \multirow{2}{*}{10} & 30 & **0.717** & 0.697 & 0.685 & 0.704 & 0.637 & 0.649 & 0.599 \\  & & 50 & **0.686** & 0.660 & 0.670 & 0.668 & 0.642 & 0.645 & 0.584 \\  & & 75 & **0.672** & 0.611 & 0.649 & 0.655 & 0.638 & 0.598 & 0.590 \\ \cline{2-10}  & \multirow{2}{*}{20} & 30 & **0.713** & 0.720 & 0.690 & 0.700 & 0.715 & 0.680 & 0.642 \\  & & 50 & **0.698** & 0.661 & 0.718 & 0.685 & 0.636 & 0.607 & 0.615 \\  & & 75 & **0.685** & 0.681 & 0.659 & 0.646 & 0.638 & 0.638 & 0.616 \\ \hline \multirow{2}{*}{250} & \multirow{2}{*}{10} & 30 & **0.720** & 0.677 & 0.703 & 0.700 & 0.667 & 0.706 & 0.630 \\  & & 50 & **0.705** & 0.675 & 0.669 & 0.676 & 0.626 & 0.646 & 0.597 \\  & & 75 & **0.709** & 0.670 & 0.680 & 0.663 & 0.617 & 0.622 & 0.601 \\ \cline{2-10}  & \multirow{2}{*}{20} & 30 & **0.727** & 0.684 & 0.709 & 0.705 & 0.660 & 0.655 & 0.628 \\  & & 50 & **0.718** & 0.683 & 0.675 & 0.682 & 0.690 & 0.708 & 0.668 \\  & & 75 & **0.710** & 0.677 & 0.686 & 0.669 & 0.658 & 0.648 & 0.618 \\ \hline \hline \end{tabular}
\end{table}
Table 2: The results of classification accuracy over 50 repeated experiments in Setting 2.
on using Setting 2 as a data generation process in order to mimic the scenario with temporal features in neural activity experiments.

In the node-wise interpretation tasks, our goal is to recover as many important nodes in the synthetic data as possible. This is particularly crucial in practice for achieving a parsimonious model that simultaneously maintains interpretability and performance. The ultimate aim is to exclude all non-important nodes and retain all important nodes after node-wise reasoning.

In the edge-wise interpretation tasks, we initially define the notions of important and redundant edges (REs). An important edge refers to an edge connecting two important nodes. In contrast, all other edges not satisfying this condition are considered non-important edges. In the task, we seek to minimize the redundant edges in the trained MaGNet estimation model. To evaluate the interpretation model performance, we define two metrics: the absolute metric (AM) and the relative metric (RM), as follows:

\[\text{AM}:=\frac{\text{\# of existing RE after reasoning}}{\text{\# of all possible RE}}, \tag{8}\]

and

\[\text{RM}:=\frac{\text{\# of existing RE before reasoning}-\text{\# of existing RE after reasoning}}{\text{\# of existing RE before reasoning}}. \tag{9}\]

In the feature-wise interpretation tasks, we have different interpretations over different settings. In Setting 1, the interpretation task is identical to the variable selection on the node features, as we aim to identify the significant dimensions of the feature contributing to the classification. In contrast, in Setting 2 where the node feature is temporally generated, the task aims to detect an important time window as in the rat hippocampus experiments.

#### 6.2.2 Results on Interpretation Tasks

In this section, we present the results of the three types of interpretation task experiments. For the purposes of comparison, we consider two competing methods: the modified GroupLasso (Meier et al., 2008) and IntGRAD (Sundararajan et al., 2017), a gradient-based model explanation approach designed for deep neural networks. To implement GroupLasso for the node-wise and edge-wise interpretation tasks, we regard each node as a variable in a group and perform
variable selection for the important nodes. For the implementation of IntGRAD, we adapt the module of IntGRAD to the graph neural network settings and apply it to our trained MaGNet estimation model. Note that for the feature-wise interpretation tasks, we only evaluate the model performance of the MaGNet interpretation method, because the two competing methods are not able to perform such type of task.

The results of the node-wise interpretation tasks are presented in Table 3 over 50 repeated experiments. Our results reveal several key findings. First, the recovery rate of important nodes using the MaGNet interpretation model consistently outperforms the competing methods across diverse settings. These advantages are due to the following strength of our method: it leverages the information gain to directly assess the reduction of uncertainty on the node subgraph. This technique incorporates statistical uncertainty as a measurement criterion instead of applying either fully deterministic gradient-based methods or variable selection methods. Another advantage of our method is the reparametrization strategy for continuously approximating discrete variables, which makes the proposed interpretation framework more computationally stable than the

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**Metric** & **Important Nodes** & **All Nodes** & **MaGNet** & **GroupLasso** & **IntGRAD** \\ \hline AM & 10 & 30 & **0.162** & 0.268 & 0.234 \\  & & 50 & **0.187** & 0.309 & 0.267 \\ \cline{2-6}  & 20 & 30 & **0.132** & 0.245 & 0.194 \\  & & 50 & **0.157** & 0.273 & 0.222 \\ \hline RM & 10 & 30 & **0.811** & 0.741 & 0.764 \\  & & 50 & **0.783** & 0.702 & 0.746 \\ \cline{2-6}  & 20 & 30 & **0.854** & 0.765 & 0.803 \\  & & 50 & **0.809** & 0.731 & 0.766 \\ \hline \hline \end{tabular}
\end{table}
Table 4: The edges-interpretation performance over two metrics with 50 repeated experiments, where lower AM is better, and higher RM is better.

\begin{table}
\begin{tabular}{c c c c c c} \hline \hline
**Sample Size** & **Important Nodes** & **All Nodes** & **MaGNet** & **GroupLasso** & **IntGRAD** \\ \hline
100 & 10 & 30 & **0.874** & 0.720 & 0.801 \\  & & 50 & **0.845** & 0.684 & 0.736 \\ \cline{2-6}  & 20 & 30 & **0.882** & 0.738 & 0.811 \\  & & 50 & **0.868** & 0.704 & 0.773 \\ \hline
250 & 10 & 30 & **0.887** & 0.729 & 0.813 \\  & & 50 & **0.856** & 0.692 & 0.746 \\ \cline{2-6}  & 20 & 30 & **0.894** & 0.748 & 0.821 \\  & & 50 & **0.877** & 0.715 & 0.785 \\ \hline \hline \end{tabular}
\end{table}
Table 3: The node-wise interpretation performance over 50 repeated experiments.
competing methods. Also, due to the fact that the MaGNet interpretation model is particularly designed for the graph neural network, it inherits the properties of the trained MaGNet estimation model, which is able to leverage local and global information to do the interpretation.

The results of the edge-wise interpretation tasks are summarized in Table 4. It shows that our proposed method has a high edge reduction rate, underlining its effectiveness in pruning the graph of non-important edges. Importantly, the method exhibits the ability to retain significant edges, suggesting an inherent proficiency in discriminating between important and non-important edges and, thus, preserving the true graph structure. This performance is consistently validated across different settings, demonstrating the proposed method's robustness and adaptability to varying graph sizes. In addition, our model identifies and removes redundant edges while maintaining critical connections between significant nodes. This results in a more parsimonious and interpretable graph structure in practice.

The results of the feature-wise interpretation task are reported in Figure 2. It indicates that our model achieves a high rate of identifying key features or time windows, and successfully detecting the most influential time windows across various settings.

Figure 2: The identified important time window with different temporal feature dimensions.
Application to Local Field Potential Activity Data from the Rat Brain

In this section, we apply the proposed method to brain activity data recorded from an array of electrodes implanted inside the brain. The brain region of interest is the hippocampus, a region near the middle of the rat brain known to be important for the temporal organization of our memories and behaviors. Although it is well established that the hippocampus plays a key role in this function across mammals, the underlying neuronal mechanisms remain unclear. To shed light on these underlying mechanisms, we previously recorded neural activity in the hippocampus of rats performing a complex sequence memory task (Allen et al., 2016) (as such high-precision data are currently not available in humans). Using that dataset, our objective here is to apply the proposed method to identify key functional relationships in the local field potential (LFP) activity simultaneously recorded across electrodes during task performance, as this information could provide novel insights into potential functional networks within that region.

The LFP activity data were collected from the CA1 region of the hippocampus while rats performed an odor sequence memory task (Fig 3). In this task, rats received repeated presentations of odor sequences (e.g., ABCDE) at a single odor port and were required to identify each item as either "in sequence" (InSeq; e.g., ABC...) or "out of sequence" (OutSeq; e.g., ABD...). Importantly, the recordings were performed from surgically implanted electrodes (tetrodes), organized into two bundles, which spanned much of the proximo-distal axis of dorsal CA1. This experimental design thus provides a unique opportunity to directly examine the anatomical distribution of information processing along that axis.

In recent work, Shahbaba et al. (2022) showed that information about trial content, such as the identity of the odor presented and whether it was presented in or out of sequence, could be accurately decoded from the ensemble spiking activity. However, that study did not determine whether task-relevant information was also contained in the local field potential activity. A fundamentally different data type from the discrete neural spiking activity, the LFP's continuous signal is more challenging to decode. To our knowledge, there are only two reports which exclusively use LFP to successfully decode spatial information in the hippocampus, which require
high-density recordings (Taxidis et al., 2015; Agarwal et al., 2014), and none showing decoding of nonspatial information from hippocampal LFP alone. To address this gap in knowledge, here we examined whether the content of odor trials can be decoded from hippocampal LFP activity and, if so, whether the pattern varies over space (electrodes) and time.

For this analysis, we focused on decoding the two main trial types (InSeq and OutSeq) using LFP activity from the 0-500 ms period (0 = odor onset), a time period in which there are no overt differences in the behavior of the animals between InSeq and OutSeq trials. We considered each rat's data an independent dataset and performed the classification evaluation task separately. For each rat's data, we randomly selected 200 graph instances as the training set and the other 30 graph samples as the testing set. Figure 4 shows that the MaGNet estimation model achieves

Figure 3: **(a.)** The task involves repeated presentations of sequences of odors and requires rats to determine whether each odor was presented “in sequence” (InSeq; e.g., ABC\(\dots\)) or “out of sequence” (OutSeq; e.g., AB\(\underline{\text{D}}\dots\)). Using an automated delivery system (left), all odors were presented in the same odor port (median interval between odors \(\sim\)5 s). Recordings were performed from electrodes organized into two bundles (right), which spanned much of the proximo-distal axis of dorsal CA1. **(b.)** In each session, the same sequence was presented multiple times, with approximately half the presentations including all InSeq trials (left) and the other half including one OutSeq trial (right). Each odor presentation was initiated by a nosepoke and rats were required to correctly identify each odor as either InSeq (by holding their nosepoke response until a tone signaled the end of the odor at 1.2 s) or OutSeq (by withdrawing their nose before the signal; <1.2 s) to receive a water reward. Incorrect responses resulted in termination of the sequence. **(c.)** Location of three electrode tips (red circles). The leftmost and rightmost electrodes approximate the extent of the CA1 transverse axis recorded in each animal
the best performance for all the rats. In particular, the proposed method indicates improvements of 8.4%, 6.9%, and 8.9% when compared to existing graph neural networks, such as GAT, GIN, and GCN, respectively. The improvements are due to our method's integration of both low-order and high-order information, effectively utilizing more information into the latent representation. In addition, by removing the non-linear feature transformation and the self-loop augmentation in GNNs, our method is less likely to suffer the _over-smoothing_ and _memoryless_ issue. Compared to the baseline classification methods, the improvement of the proposed method is more than 16.9% for DeepNet, which is the best among all the baseline approaches.

We then investigated the temporal dynamics of this decoding during trial periods by applying our MaGNet interpretation model. Specifically, we examined the most informative time bins for the InSeq/OutSeq classification in the first 500 ms of trials in Figure 4(a). We found that most significant time bins occurred between \(\sim\)180 ms and \(\sim\)320 ms after the rats poked into the port. This timeline is consistent with reports of hippocampal neurons responding to odor information in as little as 100ms (Allen et al., 2020) and with the expected timeline of InSeq/OutSeq identification within trials. This implies that the MaGNet interpretation model successfully identifies the important time window.

Figure 4: Barplot of estimation accuracy for the MaGNet estimation model and alternative competing approaches on decoding the two main trial types.
In addition, we found that most informative electrodes clustered in the distal region of CA1 in Figure 5b. In fact, across all 5 animals, the majority (86.7%) of significant electrodes were in distal CA1, and more than half of all electrodes in distal CA1 reached significance. We further

Figure 5: **(a.)** Significant decoding of InSeq and OutSeq trials based on LFP activity during the first 500ms of odor trials. Scores peak during the 185-320 ms period, prior to the behavioral response. Grey traces indicate individual subject decodings, the black line indicates the mean across subjects. **(b.)** Informative electrode nodes and edges clustered in the distal region of CA1. Schematic showing side view of electrode bundles implanted across the CA1 proximal-distal axis (Top). Schematic showing a top view of the anatomical distribution of electrodes across subjects based on electrode tract reconstruction (bottom). Yellow indicates significant nodes (electrodes). Edges indicate significant relationships between electrodes. **(c.)** The clustering of informative nodes in distal CA1 is consistent with known anatomical differences in input connections. Odor information enters the hippocampus primarily through the LEC (lateral entorhinal cortex; magenta), which more strongly projects to the distal segment of CA1. In contrast, the MEC (medial entorhinal cortex; blue) more strongly projects to proximal CA1. Approximate location of the implanted electrode bundles are shown.
found that the majority of significant edges found by the MaGNet interpretation model also clustered in distal CA1. This distribution of significant nodes and clusters suggests distal CA1 plays a more important role in representing InSeq/OutSeq information than proximal CA1, a pattern consistent with known differences in their anatomical connections (Figure 5c). However, the observation that a number of significant edges also extended into proximal CA1 suggests that functional interactions among the two segments of CA1 are critical for task performance.

Furthermore, the distal region of the hippocampus involves a more explicit layer-wise division. Specifically, it is composed of three neuronal function segments: Stratum Pyramidale (SP), Stratum Radiatum (SR), and Stratum Oriens (SO). These segments play a crucial role in the intricate circuitry of the hippocampus, encompassing not only the pyramidal neurons but also various types of interneurons that regulate the activity of these pyramidal neurons. The connectivity and interactions among neurons in these layers support complex neural conduction that contributes to the formation and retrieval of memories. As an illustrating example, Figure 6 showcases that the influential nodes identified by the MaGNet interpretation model are consistently located across all three segments, suggesting a pattern of co-activation and cross-functional firing among the neurons in these segments, essential for processing the trisynaptic circuit (Tsao et al., 2018). This co-activation and cross-firing pattern also affirms the capability of our estimation model to incorporate high-order information.

Figure 6: Three neuronal function segments in the distal area of the hippocampus. The green, red, and orange shadow area indicates the segment Stratum Pyramidale, Stratum Radiatum and Stratum Oriens, respectively.
Discussion

In this paper, we have proposed a novel graph neural network framework, MaGNet, which is able to effectively integrate the various-order information from both low-order and high-order information to allow powerful latent representation. Furthermore, MaGNet includes an interpretation component, which offers a tractable framework for identifying the influential subgraphs including important node, edge, and node features. In addition, we have also established rigorous theoretical foundations to assess the efficacy, statistical complexity, and generalizability of the MaGNet model. These theoretical results ensure that the proposed model is reliable and effective, contributing to the practical utility of MaGNet in various applications.

We applied the MaGNet framework to LFP activity data recorded from the hippocampus of rats as they performed a challenging non-spatial sequence memory task. Using this model, we were able to decode the trial type (whether the odor was presented in or out of sequence) as well as identify the most informative trial periods and electrodes. Therefore, not only did the model provide the first direct evidence of decoding non-spatial trial content from hippocampal LFP activity alone, it also provided a high degree of specificity about how this information was distributed over space and time. This neuroscience result is consistent with a growing literature on the influence of anatomical gradients on information processing within brain regions (van Strien et al., 2009; Witter et al., 2017; Knierim et al., 2014), specifically with evidence that inputs carrying non-spatial information more strongly project to distal CA1 than proximal CA1 (Haberly and Price, 1978; Agster and Burwell, 2009).

While our proposed framework maintains the sequential training approach as typical of most existing GNN architectures, the proposed training and optimization can still be costly due to the fusion step for multiple component models. Therefore, identifying ways to optimize these operations is a potential area for future research. Another potential direction for exploration is to extend the current framework to accommodate different types of tasks. For example, rather than solely focusing on graph classification, the framework can be extended for node classifications, link prediction, and beyond. By exploring these directions, we can continue to broaden the utility and effectiveness of our framework for a wide range of applications.
Furthermore, there is a significant need to extend the current framework to accommodate dynamic settings. While modeling time-varying changes and dynamic systems holds central importance in numerous real-world applications, the current MaGNet framework, along with the majority of GNN models, is primarily tailored for static graph data. These models are capable of incorporating structural information into the learning process, but they fall short in capturing the evolution of dynamic graphs. Typically, dynamics in a graph refer to node attribute modifications or edge-structure changes, including the additions and deletions of nodes or edges. As a possible expansion of the existing MaGNet framework, we will explore the incorporation of node and edge activation functions to signify and capture the presence of the nodes and edges within each timestamp. This will enable subsequent utilization of attention mechanisms such as self-attention and neighborhood attention, which have shown efficacy in foundational models (Bommasani et al., 2021), to account for historical time-evolved information from preceding timestamps.
