FLOWX: TOWARDS EXPLAINABLE GRAPH NEURAL NETWORKS VIA MESSAGE FLOWS
Year: 2021

Abstract:
We investigate the explainability of graph neural networks (GNNs) as a step towards elucidating their working mechanisms. While most current methods focus on explaining graph nodes, edges, or features, we argue that, as the inherent functional mechanism of GNNs, message flows are more natural for performing explainability. To this end, we propose a novel method here, known as FlowX, to explain GNNs by identifying important message flows. To quantify the importance of flows, we propose to follow the philosophy of Shapley values from cooperative game theory. To tackle the complexity of computing all coalitions’ marginal contributions, we propose an approximation scheme to compute Shapley-like values as initial assessments of further redistribution training. We then propose a learning algorithm to train flow scores and improve explainability. Experimental studies on both synthetic and real-world datasets demonstrate that our proposed FlowX leads to improved explainability of GNNs.

heading: 1 INTRODUCTION
text: With the advances of deep learning, graph neural networks (GNNs) are achieving promising performance on many graph tasks, including graph classification (Xu et al., 2019; Gao & Ji, 2019; Chen et al., 2020), node classification (Kipf & Welling, 2017; Veličković et al., 2018; Wu et al., 2019), and graph generation (Luo et al., 2021; You et al., 2018). Many research efforts have been made to develop advanced graph operations, such as graph message passing (Kipf & Welling, 2017; Veličković et al., 2018; Li et al., 2019), graph pooling (Yuan & Ji, 2020; Zhang et al., 2018; Ying et al., 2018), and 3D graph operations (Schütt et al., 2017; Klicpera et al., 2020). Deep graph models usually consist of many layers of these operations stacked on top of each other interspersed with nonlinear functions. The resulting deep models are usually deep and highly nonlinear and complex. While these complex systems allow for accurate modeling, their decision mechanisms are highly elusive and not human-intelligible. Given the increasing importance and demand for trustworthy and fair artificial intelligence, it is imperative to develop methods to open the black-box and explain these highly complex deep models. Driven by these needs, significant efforts have been made to investigate the explainability of deep models on images and texts. These methods are developed from different perspectives, including studying the gradients of models (Simonyan et al., 2013; Smilkov et al., 2017; Yang et al., 2019), mapping hidden features to input space (Zhou et al., 2016; Selvaraju et al., 2017), occluding different input features (Yuan et al., 2020a; Dabkowski & Gal, 2017; Chen et al., 2018), and studying the meaning of hidden layers (Yuan et al., 2019; Olah et al., 2018; Du et al., 2018), etc. In contrast, the explainability of deep graph models is still less explored. Since graph data contain limited locality information but have important structural information, it is usually not natural to directly extend image or text based methods to graphs. Recently, several techniques have been proposed to explain GNNs, such as XGNN (Yuan et al., 2020b), GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), and SubgraphX (Yuan et al., 2021), etc. These methods mainly focus on explaining graph nodes, edges, features, or subgraphs.
In this work, we observe and argue that message flows are the inherent functional mechanism of GNNs and thus are more natural and intuitive for studying the explainability of GNNs. To this end, we propose a message flow based explanation method, known as FlowX, to explain GNNs. FlowX attributes GNN predictions to message flows and studies the importance of different message flows. We first develop a systematic framework that lays the foundation on message flows for naturally explaining the message passing in GNNs. With our framework, the FlowX first quantifies
the importance of flows by following the phylosophy of Shapley values. Since message flows cannot be directly quantified to calculate marginal contributions, we propose an approximation scheme as the initial assessments of different flows. We then propose a learning-based algorithm taking advantages from the initial assessments to capture important flows in predictions. We conduct extensive experiments on both synthetic and real-world datasets. Experimental results show that our proposed FlowX outperforms existing methods significantly and consistently. Both quantitative and qualitative studies demonstrate that our proposed FlowX leads to improved explainability of GNNs.

heading: 2 RELATED WORK

heading: 2.1 GRAPH NEURAL NETWORKS
text: With the advances of deep learning, several graph neural network approaches have been proposed to solve graph tasks, including graph convolutional networks (GCNs) (Kipf & Welling, 2017), graph attention networks (GATs) (Veličković et al., 2018), and graph isomorphism networks (GINs) (Xu et al., 2019), etc. They generally follow a message-passing framework to learn graph node features. Specifically, the new features of a target node are learned by aggregating messages flows passed from its neighboring nodes. Without loss of generality, we consider the input graph as a directed graph with n nodes and m edges. The graph is denoted as G = (V,E), where V = {v1, . . . , vn} denotes nodes, and E = {eij} represents edges in which eij is the directed edge vi → vj . Then it can be represented by a feature matrix X ∈ Rd×n and an adjacency matrix A ∈ Rn×n. Each node vi is associated with a d-dimensional feature vector xi corresponding to the i-th column of X . The element aij in A represents the weight of eij , and aij = 0 indicates eij does not exist. For the t-th layer in GNNs, the message aggregation procedures can be mathematically written as a two-step computation as
Aggregate: St = Xt−1Ât, (1)
Combine: Xt = M t(St), (2)
where Xt ∈ Rdt×n denotes the node feature matrix computed by the t-th GNN layer and X0 = X . Here M t(·) denotes the node feature transformation function at layer t and Ât is the connectivity matrix at layer t. Note that we name the elements in Ât as layer edges and âtjk indicates the layer edge connecting node j and k in layer t. For example, in GCNs, the transformations are defined as σ(W tSt) and Ât = D− 1 2 (A+ I)D− 1 2 where W t ∈ Rdt×dt−1 is a trainable weight matrix, σ(·) denotes the activation function, I is an identity matrix to add self-loops to the adjacency matrix, and D denotes the diagonal node degree matrix. We can stack T GNN layers on top of each other to form a T -layer network, and the network function can be expressed as
f(G) = g(MT (MT−1(· · ·M1(X0Â1) · · · )ÂT−1)ÂT ).
When f(G) is a graph classification model, g(·) generally consists of a readout function, such as global mean pooling, and a multi-layer perceptron (MLP) graph classifier. Meanwhile, when f(G) is a node classification model, g(·) represents a MLP node classifier.

heading: 2.2 EXPLAINABILITY OF GRAPH NEURAL NETWORKS
text: A major limitation of GNNs is their lack of explainability. Thus, different methods have been proposed to explain the predictions of GNNs, such as GraphLime (Huang et al., 2020), GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), PGMExplainer (Vu & Thai, 2020), SubgraphX (Yuan et al., 2021), XGNN (Yuan et al., 2020b), and GraphSVX (Duval & Malliaros, 2021). These methods can be mainly grouped into four categories based on the views of their explanations. First, several techniques provide explanations by identifying important nodes in the input graph. For example, GradCAM (Pope et al., 2019) measures node importance by combining the hidden features and gradients; LRP (Baldassarre & Azizpour, 2019) and Excitation BP (Pope et al., 2019) decompose the predictions into several terms and assign these terms to different nodes; PGM-Explainer (Vu & Thai, 2020) builds a probabilistic graphical model by randomly perturbing the node features and employs an interpretable Bayesian network to generate explanations. Second, several existing methods, such as GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), and GraphMask (Schlichtkrull
et al., 2021), explain GNNs by studying the importance of different graph edges. These methods follow a similar high-level idea that learns masks to identify important edges while maximizing the mutual information. Next, the recent study SubgraphX (Yuan et al., 2021) proposes to explain GNNs via subgraphs. It incorporates the Monte Carlo tree search algorithm to explore subgraphs and employs Shapley values to measure the importance. Finally, XGNN (Yuan et al., 2020b) focuses on model-level explanations, which can provide high-level insights and general understanding. It proposes to generate graph patterns that can maximize a certain model prediction. While those methods explain GNNs from different views, none of them can provide explanations in terms of message flows and we believe message flows are more natural for performing explainability.
To the best of our knowledge, GNN-LRP (Schnake et al., 2020) is the only algorithm that explains GNNs by relevant walks. The relevant walk is defined as a T -length ordered edge sequence that corresponds to a T -step directed path on the input graph. To study walk explainability, GNN-LRP considers the GNN prediction as a function and decomposes it using higher-order Taylor expansions to distribute prediction scores to relevant walks. Specifically, by using T -order Taylor expansion with a proper root, each term in the Taylor expansion corresponds to a relevant walk and is regarded as the importance score. While our proposed FlowX shares a similar explanation target, i.e., flow/walk, with GNN-LRP, our method is fundamentally different. The GNN-LRP is developed based on score decomposition while our method follows the phylosopy of Shapley values from cooperative game theory as initial assessments and proposes a learning-based algorithm for the score generation. In addition, GNN-LRP has several constraints on the activation function and bias term used in GNNs while our method can be applied to general GNN models. Furthermore, as the GNN-LRP follows the Gradient× Input scheme, it may not pass the model parameter randomization test and may not be sensitive to model parameters (Adebayo et al., 2018).
Differences with Other Methods Using Shapley Values: Shapley values are commonly used in explaining machine learning methods. In particular, a recent study proposes a surrogate method, known as GraphSVX (Duval & Malliaros, 2021), to explain GNNs with both node and feature masks. Another recent study proposes SubgraphX (Yuan et al., 2021), which employs a search algorithm to explore and identify subgraphs with high Shapley scores. While these methods use Shapley values, there are several fundamental differences. First, our proposed method focuses on explaining message flows, which are the most basic and natural units for explanations as GNNs are based on message passing schemes. Second, we only use Shapley-like values as initial approximations to facilitating further training. Experiments show that the learning step is very important. Furthermore, the fundamental difference with GraphSVX is reflected in the fact that, while our method is a perturbation-based method, GraphSVX is a surrogate method (Yuan et al., 2020c). Due to these differences, we show in experiments that our FlowX achieves more natural and improved performance as compared with other methods.

heading: 3 THE PROPOSED FLOWX
text: While existing methods mainly focus on explaining GNNs with graph nodes, edges, or subgraphs, we propose to study the explainability of GNNs from the view of message flows. We argue that message flows are the fundamental building blocks of GNNs and it is natural to study their contributions towards GNN predictions. With our message flow framework, we propose a novel method, known as FlowX, to investigate the importance of different message flows. Specifically, we follow the phylosophy of Shapley values (Kuhn & Tucker, 1953) from game theory and propose an marginal contribution approximation scheme for them. In addition, a learning-based algorithm is proposed to improve the explainability of message flows.

heading: 3.1 A MESSAGE FLOW VIEW OF GNNS
text: We consider a deep graph model with T GNN layers. Each GNN layer aggregates 1-hop neighboring information to learn new node embeddings. Hence, for any node, the outgoing messages are transmitted within its T -hop neighbors. Then the outputs of GNNs can be regarded as a function of such transmitted T -step messages, which are named as message flows in this work. Formally, we introduce the concept of message flows and message carriers as follow:
Definition 1: Message Carrier. We use the connectivity matrix to represent the carriers for message flows. Given the connectivity matrix Ât at layer t, the layer edge âtij represents the message carrier with which the message passes from node vi to vj at layer t.
Note that we use the superscript t to distinguish the message carriers in different layers since their corresponding message flows are different. Then the set of all message carriers, i.e., all layer edges, is defined as A = {· · · , â1uv, · · · , âtuv, · · · , âTuv, · · · } and |A| = |E| × T . Definition 2: Message Flow. In a T -layer GNN model, we use Fijk...`m to denote the message flow that starts from node vi in the input layer, and sequentially passes the message to node vj , vk, . . . , v` until to node vm in the final layer T . The corresponding message carriers can be represented as {â1ij , â2jk, . . . , âT`m}.
In a T -layer GNN model, all message flows start from the input layer and end with the final layer so that their lengths equal to T . For the ease of notations, we introduce the wildcard ∗ to represent any valid node sequence and F to denote the message flow set. For example, we can use Fij∗, F∗`m, and Fij∗`m to denote the message flow sets that share the same message carrier â1ij , âT`m, or both of them, respectively. In addition, we employ another wildcard ? to denote any single node and ?{t} to represent any valid node sequence with t nodes. For example, F?{3} means the set of valid 2-step message flow with 3 nodes. Note that the following property of message flow sets also holds:
Fij∗`m = Fij∗ ∩F∗`m. (3)
The final embeddings of node vm are determined by all incoming message flows to node vm, which can be denoted as F∗m. Since the output of the GNN model is obtained based on the final node embeddings, then it is reasonable to treat the GNN output as the combinations of different message flows. Hence, it is natural to demystify GNN models by studying the importance of different message flows towards GNN predictions.

heading: 3.2 SAMPLING MARGINAL CONTRIBUTIONS AS INITIAL ASSESSMENTS
text: While explaining GNNs with message flows seems to be promising, it is still crucial to properly measure the importance of those message flows. Hence, in this work, our FlowX proposes to follow the philosophy of Shapley value (Kuhn & Tucker, 1953) to use the marginal contributions in different flow sets as the initial assessments of flow importance. It is noteworthy that our initial assessments are not computing real Shapley values. Shapley value is a solution concept in cooperative game theory and used to fairly assign the game gain to different players. When considering marginal contributions in GNN explanation tasks, we treat the message flows as different players and the GNN prediction score as the total game gain. Formally, given the trained GNN model f(·) and the input graph G, we use F∗ to denote the set of all valid flows, i.e., all players in the game. Then given any flow Fk, we mathematically define the contribution of it as
φ(Fk) = ∑
P⊆F∗\{Fk}
W (|P |)(f ( P ∪ {Fk} ) − f(P )), (4)
where W (|P |) is a weight function assigned to each term in the summation; P denotes the possible coalition group of players, and φ(·) denotes the flow score. Here f ( P ∪ {Fk} ) − f(P ) is the
marginal contribution of flow Fk for a particular coalition group, which can be computed by the prediction difference between combining Fk with the coalition group P and only using P . Note that Eq. 4 is equivalent to the classic Shapley value when W (|P |) = |P |!(|F∗|−|P |−1)!|F∗|! where | · | is the set size. To compute the flow score φ(Fk), we need to enumerate all possible coalition groups and considers different interactions among players. However, it is time-consuming to consider all possible coalition; thus, we sample several marginal contributions to approximate the final flow score. Note that it is not possible to compute importance score with Eq. (4) for message flows since we cannot remove individual message flows from the GNN model and the finest component we can directly remove is the layer edge. For example, when removing the layer edge â1ij , the whole flow set Fij∗ is removed from the model. In addition, enumerating all possible coalition groups is time-consuming when the input graph is large-scale and the GNN model is deep. Hence, in our FlowX, we propose an approximation scheme to compute Eq. 4 based on Monte Carlo (MC) sampling (Štrumbelj & Kononenko, 2014).
Algorithm 1 INITIAL APPROXIMATIONS OF FLOW IMPORTANCE SCORES. 1: Given a trained GNN model f(·) and an input graph, the set of all layer edges is represented
as A. For each message flow Fk, two |A|-dimensional vectors S(Fk) and C(Fk) denote its importance scores and removing index counts respectively. In addition, Â denotes the set of removed layer edges and F̂ is the set of removed message flows. 2: Initialize S(·) and C(·) as zeros for all message flows. 3: for step i from 1 to Monte Carlo sampling step M do 4: Initialize the removed sets as empty that Â = ∅ and F̂ = ∅. 5: Randomly shuffle and permute the layer edge set A, denoted as Aπ(i). 6: for j from 1 to |A| do 7: Select the j-th layer edge in Aπ(i), denoted as ât`m. 8: Block the layer edge ât`m in GNN model f(·), then the removed flows are F̂ j =
F?{t−1}`m?{T−t} \ (F̂ ∩F?{t−1}`m?{T−t}). 9: Compute the prediction difference that sj = f(A \ Â)− f(A \ (Â ∪ ât`m)).
10: Update Â = Â ∪ ât`m and F̂ = F̂ ∪F?{t−1}`m?{T−t}. 11: Compute averaged score that s̄j = sj/|F̂ j |. 12: for each flow Fk in F̂ j do 13: Update S(Fk)[j] = S(Fk)[j] + s̄j . 14: Update C(Fk) = C(Fk) + 1 15: end for 16: end for 17: end for 18: For each message flow Fk, compute our marginal contribution vector S(Fk) = S(Fk)/C(Fk)
Formally, let M denote the total number of MC sampling steps. For the i-th sampling step, we use π(i) to represent a random permutation of |A| = |E| × T elements. Then the set of all layer edge A is permuted based on π(i), denoted as Aπ(i) that Aπ(i)[j] = A[π(i)[j]]. In each step i, we iteratively remove one layer edge from Aπ(i) following the order, and compute the marginal contribution. An illustration of our proposed sampling algorithm is shown in Figure 1. We use Â to denote the set of removed layer edges and F̂ to denote the set of removed message flows, which are initialized as Â = ∅ and F̂ = ∅ in the beginning of each MC sampling step. For the iteration j in sampling step i, the j-th element of Aπ(i) is removed and we assume the removed layer edge is ât`m. Then the computation operations can be mathematically written as
sj = f(A \ Â)− f(A \ (Â ∪ ât`m)), (5) Â = Â ∪ {ât`m}, (6) F̂ j = F?{t−1}`m?{T−t} \ (F̂ ∩F?{t−1}`m?{T−t}), (7)
F̂ = F̂ ∪F?{t−1}`m?{T−t}, (8)
where F̂ j denotes the removed message flows in iteration j by removing ât`m. Note that it is not equivalent to F?{t−1}`m?{T−t} since the flows in F?{t−1}`m?{T−t} may be already removed in the previous iterations. Then the score is averaged that s̄j = sj/|F̂ j | and we assign s̄j to each flow in F̂ j . By repeating such operations until all layer edges are removed, we can obtain |A| marginal scores and assign them to the corresponding flows. Note that the order information in Aπ(i) is important since in the earlier iterations, the removed flows F̂ j are interacting with a larger coalition group A \ (Â ∪ ât`m) and in the later iterations the coalition groups are small. Altogether, the steps of our marginal contribution sampling are shown in Algorithm 1.
After M MC sampling steps, our method explores M permutations of layer edge set and each flow is sampled to obtain M marginal contributions. For each flow Fk, we use a |A|-dimensional vector to store its marginal contributions, denoted as S(Fk), where S(Fk)[j] is the importance score obtained when Fk is removed in the iteration j. If the flow is removed mutiple times in the iteration j for different MC sampling steps, the scores are averaged, which also means scores obtained from coalitions with same sizes are weighed equally. Then it raises the question that how to convert S(Fk) to the final importance score of Fk. The final importance score will be a Shapley value approximation if we simply compute the summation of the elements in S(Fk). While Shapley values are promising, the approximations may not be accurate since M is always set as M << |A|! for the sake of computation efficiency. Moreover, the importance of the early iterations when more layer edges exist in the model should be different to the importance of later iterations when only few layer edges exist in diverse situations. Hence, directly computing the summation of the scores in S(Fk) may not be suitable.

heading: 3.3 LEARNING IMPORTANCE SCORES
text: With our obtained score vectors S(Fk), we propose to regard these values as initial assessments of flow importance and learn an associated importance scores with designed redistribution trick. We target to learn weights for combining the elements in S(Fk) to obtain the final importance score for each flow, which means to learn W (P ) in Eq. 4.
Formally, given the input graph G and GNN model f(·), we obtain marginal contribution vector S(Fk) for each message flow. Then the only trainable weight vector w is initialized to all equal to 0.5 with random noise, which shares the same dimension of |S(Fk)|. Note that the weights w are shared among all message flows in the same graph. Then S(Fk) can be combined via a dot product with w. Next, the flow importance is converted to layer edge importance by simply summing up the scores of all flows sharing a particular layer edge as the message carrier. Mathematically, it can be written as
s(Fk) = S(Fk) · w ∀Fk ∈ F∗, (9) s(âtuv) = ∑ Fj∈F̃ s(F j) (10)
where s(·) denotes the importance score for layer edges and flows; F̃ = F?{t−1}uv?{T−t}. Then based on s(âtuv), we can obtain a mask indicating the importance of different layer edges, denoted as M:
M = g(ŝ), (11)
where ŝ denotes all layer edges’ importance scores. g(·) is defined as an exponential redistribution operation that including an input normalization, an element-wise exponential scaling ge(x) = xr (r is a hyper-parameter), and an output normalization. We provide more intuitions in Appendix A.3. Note that the mask is normalized so that each element is in the range [0, 1]. By applying the mask to layer edges, important layer edges are restricted and the model prediction becomes
ŷ = f(Combine(G, 1−M)), (12)
where ŷ is the prediction vector and Combine(G, 1−M) means layer edges are masked out based on the values ofM. Intuitively, if important layer edges are restricted, then the prediction should change significantly. Hence, ŷ is encouraged to be different from the original prediction by learning proper
weights w. Specifically, we employ the negative cross-entropy loss that can be formally expressed as
L(ŷ, y) = m∑ c=1 1{y = c} log ŷc, (13)
where y is a scalar representing the predicted class of the original graph G, m is the number of classes, 1{·} denotes the indicator function, and ŷc is the new predicted probability for class c. After training, the final importance scores for different flows can be obtained via Eq. (9).

heading: 4 EXPERIMENTAL STUDIES

heading: 4.1 DATASETS AND EXPERIMENTAL SETTINGS
text: Datasets. We employ seven different datasets to demonstrate the effectiveness of our proposed FlowX with both quantitative studies and qualitative visualization results. These datasets are BAShapes (Ying et al., 2019), BA-LRP (Schnake et al., 2020), ClinTox (Wu et al., 2018), Tox21 (Wu et al., 2018), BBBP (Wu et al., 2018), BACE (Wu et al., 2018), BA-INFE, and Graph-SST2 (Yuan et al., 2020c), which include both synthetic and real-world data. First, BA-Shapes is a node-classification synthetic dataset that is built by attaching house-like motifs to the base Barabási-Albert graph where the node labels are determined by their own identifies and localizations in motifs. Then, BA-LRP is a graph-classification synthetic dataset that includes Barabási-Albert graphs and the two classes are node-degree concentrated graph and evenly graph. Next, ClinTox, Tox21, BBBP, and BACE are real-world molecule datasets for graph classification. The chemical molecule graphs in these datasets are labeled according to their chemical properties, such as whether the molecule can penetrate a blood-brain barrier. Finally, Graph-SST2 is a natural language sentimental analysis dataset that converts text data to graphs. These graphs are labeled by their sentiment meanings. Note that the BA-INFE dataset is a synthetic dataset for Accuracy comparisons and the results are reported in the appendix. We also report the properties and statistics of these datasets in Appendix B Table 2.
GNN Models. In our experiments, we consider GCNs (Kipf & Welling, 2017) and GINs (Xu et al., 2019) as our graph models for all datasets. We adopt 2-layer GNNs for node classification and 3-layer GNNs for graph classification. The graph models are trained to achieve competitive performance and the details are reported in Table 2.
Baselines. With the trained graph models, we quantitatively and qualitatively compare our FlowX with eight baselines, including GradCAM (Pope et al., 2019), DeepLIFT (Shrikumar et al., 2017), GNNExplainer (Ying et al., 2019), PGExplainer (Luo et al., 2020), PGMExplainer (Vu & Thai, 2020), SubgraphX (Yuan et al., 2021), GNN-GI (Schnake et al., 2020), GNN-LRP (Schnake et al., 2020). Note that since these methods cannot be directly compared, we set the explanation target to graph edges for fair comparisons that the explanations are converted to edge importance scores if needed. More implementation details about explanation methods setting and GNN models can be found in the supplementary material. We used the datasets and implementations of the comparing algorithms in the DIG library (Liu et al., 2021). We will release our code after the anonymous review period.

heading: 4.2 QUANTITATIVE STUDIES
text: We first quantitatively compare different explanation methods. We follow the existing studies (Yuan et al., 2020c; Pope et al., 2019) and employ two metrics to evaluate the explanations: Fidelity and Sparsity. Good explanations should be faithful to the model and capture the discriminative features for the predictions. When such input features are removed, the original predictions should change significantly. The Fidelity score measures the change of predicted probabilities when removing important input features identified by different explanation methods. Higher Fidelity scores indicate the removed features are more important to the predictions and hence the explanations are more faithful to the model. In addition, the other desired property of explanations is sparsity. To encourage the explanations to be more human-intelligible, they should contain fewer but more important features. Hence, we also employ the Sparsity metric which measures the percentage of input features that are identified as important. Higher Sparsity scores indicate that fewer features are identified as important in the explanations. The formulations and details of these metrics are discussed in the appendix C.
It is noteworthy that these two metrics are highly correlated since the predictions tend to change more when more input features are removed. Then explanations with higher Sparsity scores tend to have lower Fidelity scores. Hence, we argue that the Fidelity scores need to be compared under a similar Sparsity level. Specifically, we control the Sparsity scores of explanations and compare the corresponding Fidelity scores. For each dataset, we randomly select 100 samples and conduct such quantitative evaluations. The results are reported in Figure 2 where we show the plots of Fidelity scores with respect to different Sparsity levels. Obviously, FlowX performs better on all of the real-world datasets that FlowX consistently achieves higher Fidelity across all Sparsity levels. Meanwhile, the results on synthetic datasets are competitive when the Sparsity is low and our FlowX still performs significantly better with high Sparsity levels. Note that our method is shown to achieve stable performance on different datasets and GNN models, which is promising for the applications and generalizations of our FlowX. Meanwhile, our proposed method, GNN-LRP, and GNN-GI perform better than the other methods, which indicates the superiority of the methods based on flows and walks. Note that GNN-LRP achieves good results on synthetic datasets but surprisingly, it is not better than GNN-GI on real-world datasets. In addition, PGExplainer only obtains competitive results on simple synthetic datasets but not on complex real-world datasets. Note that, for SubgraphX, its subgraph explanations are not directly comparable so we convert its subgraph explanations to edge explanations. It performs well on the BA-LRP and BBBP datasets but not as expected on the other datasets. We believe the reason may be that such converting destroys the continuousness and completeness of its subgraph explanations. Last but not least, we wish to mention that the comparisons using edge actually limit the performance of the flow-based methods due to explanation granularity changes. More results and time comparisons are also reported in the appendix.

heading: 4.3 ABLATION STUDIES
text: To demonstrate the effectiveness of our proposed learning-based algorithm in Section 3.3, we study the performance when only using the average of our Shapley initial assessments as the importance scores without further learning, denoted as FlowX∗. In addition, to show the Shapley initial assessments are
necessary, we study the method only using the learning algorithm with randomly initialized scores for different flows, denoted as FlowX†. We compare our FlowX with these two baselines on GCN models, and the results are reported in Table 1. Obviously, for all datasets, the FlowX and FlowX∗ outperform FlowX† significantly, showing that the Shapley initial assessments are appropriate and necessary. In addition, for synthetic datasets, the FlowX and FlowX∗ show competitive performance, which indicates our initial assessments can obtain reasonable results for simple models and datasets. However, FlowX obtains significantly better performance on complex real-world datasets than FlowX∗, showing the effectiveness of our proposed learning-based algorithm. More experimental details will be attached in the supplementary material.

heading: 4.4 VISUALIZATION OF EXPLANATION RESULTS
text: Finally, we report the visualization results of different explanation methods in Figure 3. Specifically, we show the explanations of a molecule graph from ClinTox dataset. The smiles string of this input graph is C(CC(C(F)F)(C(=O)[O-])[NH3+])C[NH3+] and the model to be explained is a 3-layer GCN model. In this input graph, motifs 1 and 2 form a Carboxylate Anion; motif 4 and 5 both contain an Ammoniumyl while motif 3 corresponds to a CF2 group. It is clear that our FlowX finds all motifs in the explanation and obtains the highest Fidelity score. It indicates that those motifs are indeed important for the model to make predictions. In addition, we can conclude that the model combines these motifs for the predictions since we generally observe that explanations with more motifs chosen tend to have higher Fidelity scores. Furthermore, we find that DeepLIFT identifies all motifs but focuses more on the self-loops of atoms. However, its Fidelity score is much lower than GNNExplainer who only identifies four motifs but focuses on the information transmissions between atoms. Hence, we believe the interactions among different atoms contribute more to the model predictions, which further indicates the Shapley value is a promising solution for studying the explainability of GNN models. More visualization results are reported in the supplementary material.

heading: 5 CONCLUSIONS
text: We study the explainability of deep graph models, which are generally treated as black-boxes. From the inherent functional mechanism of GNNs, we propose FlowX to explain GNNs by studying message flows. Our FlowX first computes Shapley values approximations as initial assessment and then incorporates a learning refinement. Extensive experiments demonstrate our FlowX achieves significantly improved explanations.

heading: REPRODUCIBILITY STATEMENT
text: Several efforts are made to ensure the reproducibility of our work. To demonstrate the details of proposed algorithm, we provide the rigorous description of our approach’s first stage in Algorithm 1, the second stage learning process in Subsection 3.3 and the concrete setting in appendix B.4. For experiments, besides the detailed description of datasets and evaluation metrics in Section 4, we provide model configurations, training details, target mapping, and metrics details in appendix B and appendix C. Our implementations will be publicly available once the paper is published.

heading: A ALGORITHM ANALYSIS & INTUITIONS

heading: A.1 MESSAGE FLOW INTUITION
text: Message flows can provide multi-hop dependence association, i.e., the reason of a given node why it is labeled as susceptible will be explained by a message flow beginning from an affected node to the given node in a virus infection dataset. Therefore, the explanation is that the message flow transmits the affection information from the affected node to the given node making it a susceptible one. The upper figures of Figure. 4 show that how a node’s message can be transmitted from the node to which node across which path.

heading: A.2 COMPLEXITY ANALYSIS
text: Message flows. Considering the total number of message flows |F |, its loose upper bound is |F | = O(|E|T ). With the consideration of connectivity, given the largest outgoing degree of nodes in the graph d+, the tighter upper bound of the number of message flows is |F | = O(|E|(d+)T−1). Marginal contribution sampling. It is noticeable that the deep model’s forward operations are the most time-consuming operations, denoted as O(Tf ). Without any parallel consideration, the time complexity of marginal contribution sampling is Tmcs = O(M |A|Tf ). When we consider parallel implementations, the most simple improvement is to use the GPUs acceleration to move line 9-15 out of the two loops and to make line 8 executed in O(1), which leads to the time complexity Tmcs = O(M |A|+Tf ). The space complexity should be considered because there is a big flow score matrix Scpu = O(M |A||F |). As for the graphic memory requirement, this parallel implementation consumes Sgpu = O(M |A|Sf ), where O(Sf ) = O((|V |+ |E|)dT +Sp) is the deep model memory complexity; here, |V | = n denotes the number of nodes; Sp is the model’s parameter memory complexity. In our implementation, there is a trade-off that we only parallelly calculate the model outside the innermost loop, so that our time complexity is Tmcs = O(M(|A|+ Tf )) and the space complexity of GPUs is Sgpu = O(|A|Sf ).

heading: A.3 TRAINING DILEMMA & EXPONENTIAL REDISTRIBUTION
text: We observe that the training is challenging as it often even leads to a worse result after training. The difficulty comes from that the training target is unusual. We not only require the important layer edges’ scores to be trained to be their high-value local optimums with correct rankings but also require that all of the unimportant ones, not including those negative contribution elements, are trained to be low values because of the ranking requirement. The problem is that because the unimportant layer edges exert nearly no impact on the model results, their scores are trained to be the values that are largely random, which causes an unimportant layer edge’s score may be larger than an important one because it is random (it may be arbitrary high or low). This is not what we expected. In brief, these random scores severely interfere with the ranking of scores of the important. Therefore, exponential redistribution is designed to solve this problem as shown in Eq. 11, where the gradients of the exponential scaling ge(x) = xr decide that the high scores can be sensitively trained while the low scores are likely to be trapped in a low-value zone with gradients closed to zeros. Intuitively, considering the simplest 2-order exponential scaling ge(x) = x2 = x · x, this function can be regarded as using x as an attention map to mask itself so that high scores become easier to train while it is contrary to the low scores. It is noticeable that the use of exponential redistribution requires a relatively reasonable initialization, in order to give important elements a lower possibility to be initialized into the low-value zone, which leads to a longer training time. Hence, we initialize scores according to the marginal contribution results we obtain from the last stage so that unimportant factors are set into the low-value zone and will be likely to be trapped there having no interference to the high score ranking, while important layer edges are normally trained. Note that the final scores are not refined Shapley values, because the Shapley-like initialization only serves as heuristic values for the starting of exponential redistribution. Worth mentioning that the use of regularization terms like the size of masks and the information entropy does not work, and the use of Gumbel softmax for encouraging discreteness is useless either because both of them cannot solve the real problem.

heading: B EXPERIMENTAL SETTINGS

heading: B.1 DEEP GRAPH MODELS
text: We first introduce the details of the GNN models we try to explain. For all GNN models, we employ two message passing layers for node classification (BA-Shapes (Ying et al., 2019)) and three for graph classification tasks (ClinTox, Tox21, BBBP, BACE (Wu et al., 2018), BA-LRP (Schnake et al., 2020), Graph-SST2 (Yuan et al., 2020c)). In addition, the final classifier consists of two fully-connected layers. For graph classification tasks, the average pooling is used to convert node embeddings to graph embeddings. In addition, the feature dimensions of message passing layers and fully-connected layers are set to 300. We apply the ReLU function as the activation function after each message passing layer and fully-connected layer. In addition, the dropout is applied between two fully-connected layers. Note that we consider both GCNs and GINs for all datasets. For the GCN layer, we employ the original normalized Laplacian matrix Ât = D− 1 2 (A+ I)D− 1 2 . For the GIN layer, two fully-connected layers with the ReLU function are employed as the multilayer perceptron (MLP).

heading: B.2 TRAINING SETTINGS
text: The models on different datasets use different learning rates and decay settings. Generally, we set the learning rate to 1× 10−3 and the learning rate decay equal to 0.5 after 500 epochs. For datasets BA-LRP and Graph-SST2, the total number of epochs is 100 while we train models for 1000 epochs on the other datasets. All datasets are split into the training set (80%), validation set (10%), and testing set (10%). All experiments are conducted using one NVIDIA 2080Ti GPU.

heading: B.3 TARGETS OF EXPLANATIONS
text: Since different techniques focus on different explainable targets, such as nodes, edges, walks, flows, etc., these methods cannot be directly compared. For fair comparisons, we convert all explainable targets to graph edges. First, GNNExplainer and PGExplainer provide edge-level explanations so their results are directly used. Second, for flow-based methods and walk-based methods, including GNN-GI, GNN-LRP, and our FlowX, we convert the flow (walk) importance to edge importance by summing the total contribution of message flows (walks) that go through a particular edge. In addition, for node-based methods such as GradCAM (Pope et al., 2019) and DeepLIFT (Shrikumar et al., 2017), the node importance is mapped to edge importance that the contribution of an edge is the averaged contribution of its connected nodes. Specially, for subgraph-based method SubgraphX (Yuan et al., 2021), we pick the explainable subgraph out, then assign the edges in this subgraph instead of nodes as the explanation. Note that the absolute values of contribution are not important since the metrics focus on their relative rankings. The Figure. 4 shows how to map from message flows to edges.

heading: B.4 FLOWX SETTINGS
text: At the first stage, the MC sampling time M is set as 30. At the second stage, we randomly initialize the weight vector of our learning refinement from the uniform distribution [0, 0.1]. During the training,
we apply a learning rate equal to 0.3 without learning rate decay and train the weight vector for 500 iterations. Meanwhile, we adopt r = 8 for exponential redistribution.
FlowX∗ & FlowX†. FlowX∗ can be implemented by replacing the weight vector w in section 3.3 with a constant vector w′ that w′ = 1/|w| where |w| is the length of w. FlowX† is implemented with a similar pipeline as GNNExplainer (Ying et al., 2019). Instead of training the weights to obtain the edge mask, we use a flow mask to indicate the flow importance scores. The flow mask can be considered as the weighted flow scores in our FlowX. In addition, we do not apply exponential scaling to FlowX†. The rest of FlowX† is the same as the learning refinement part in our FlowX.

heading: C EVALUATION METRICS
text: Mathematically, a K-class classification dataset with N samples can be represented as {((Xi,Λi), yi) | i = 1, 2, . . . , N}, where (Xi,Λi) is the input graph and yi is the ground-truth label of the i-th sample. The predicted class of a GNN model is ŷi = argmaxf(Xi,Λi). Here the output of f(Xi,Λi) is a K-length vector, in which each element denote the probability of the corresponding class. Formally, we define the algorithm’s explanation as a column vector ṁi that represents a mask on explainable targets where the number of targets is |ṁi|. In this mask, the explanation method marks each selected target as a value 1, otherwise a 0. The source of the mask is a vector ċi with the same shape that stores the contribution score of each explainable target. Intuitively, we prefer to assign targets with 1s in the ṁi if the corresponding scores in ċi are relatively high. With these notations, we introduce the metrics employed for quantitative comparison.
Fidelity. Following the metrics mentioned in the graph explainability survey (Yuan et al., 2020c), we choose the metric Fidelity to measure the faithfulness of the explainability methods to the model. Mathematically, Fidelity is defined as:
Fidelity = 1
N N∑ i=1 (f(Xi,Λi)ŷi − f ṁi(Xi,Λi)ŷi), (14)
where f ṁi(Xi,Λi)ŷi denotes the output class ŷi’s probability when masking out the edges based on ṁi. Specifically, we keep the value unchanged for the explainable targets that are marked as 1 in ṁ, while setting the other explainable targets as 0. Note that ṁi = 1− ṁ is complementary mask where 1 is an all-one vector that has the same shape as ṁ. In the evaluation period, the higher Fidelity indicates the better performance of the explainability algorithm. Intuitively, Fidelity represents the predicted probability dropping after occluding the important features. Therefore, Fidelity can measure whether the chosen features are important to the predictions from the model’s perspective.
Sparsity. In order to fairly compare different techniques, we argue that controlling the Sparsity (Pope et al., 2019; Yuan et al., 2020c) of explanations is necessary. Mathematically, Sparsity is defined as:
Sparsity = 1
N N∑ i=1
( 1− ∑|ṁi| k=1 ṁki |ṁi| ) . (15)
A higher Sparsity score means fewer explainable targets are selected. Note that the selections of targets are determined by the ranking of the contribution score in ċi. It is noteworthy that these two metrics are highly correlated since the predictions tend to change more when more input features are removed. Then explanations with higher Sparsity scores tend to have lower Fidelity scores. For fair comparisons, we compare different methods with similar sparsity levels.

heading: D EXPERIMENTAL RESULTS
text: In this section, we provide additional quantitative comparisons and visualization results.

heading: D.1 QUANTITATIVE RESULTS
text: We report the quantitative results of GCN and GIN models in Table 4 and Table 5. Obviously, our method obtains the best averaged Fidelity score for 5 out of 7 GCN models. Note that the performance is competitive for synthetic dataset BA-Shapes and BA-LRP and our method performs significantly better on the complex real-world datasets. In addition, we show the plots of Fidelity scores with respect to different Sparsity levels for all datasets and GCN models in Figure 5. Clearly, the performance of our method is stable and better.
Accuracy comparisons: To further demonstrate the effectiveness of our proposed method, we conduct experiments to compare the explanatory accuracy of three methods that usually evaluated in terms of accuracy. The results are reported in Table 6. These results are obtained from a synthetic graph classification dataset BA-INFE. Each graph includes a base BA-graph and four types of motifs. Two classes of motifs are attached (2 motifs denote the one property, 2 for another; let’s denote them as class/property A and B). We first connect the same number (1 3) of each class’s motifs to the base graph. Then we attach 1 to 3 extra motifs from one of the classes (for example, class A) to the base graph; thus, the number motifs in class A is more than class B, representing the graph tends to have the specific property (property A). The ground-truths of a graph are the elements in those motifs including edges (undirected) and self-loop connections. The number of ground-truth elements that are covered by explanation edges is denoted as the hit number. The accuracy is defined as (hit number) / (total number of ground-truth elements).
We also show the time cost of different methods in Table 3. Specifically, we report the averaged computation time for 40 graphs in the ClinTox (Wu et al., 2018) dataset. It shows that our FlowX performs the best with a reasonable time cost.

heading: D.2 EXPLANATION VISUALIZATIONS
text: We provide additional examples to visually compare the explanations by six methods. The results are reported in Figure. 6, and 7. We can conclude that the message interactions between nodes are more important than self-loops. For example, in Figure. 6, all of DeepLIFT, GNNExplainer, and FlowX can identify the motifs; however, FlowX attributes more on the node interactions and obtain a higher Fidelity score.
Figure 6: Sample explanation results of different methods on Tox21.
Figure 7: Sample explanation results of different methods on ClinTox.

Review summaries:
FlowX is an explanation method for Graph Neural Networks. It derives importance measures for message flows, which are inherent function mechanisms of GNNs. These measures are initialized with a Monte Carlo approximation of the Shapley Values from Game Theory, applied to message flows. They are refined by learning a linear transformation, before being converted back to layer edge importance score ultimately provided as an explanation. Various experiments are performed to assess the relevance of the proposed explainer.

This paper explains GNNs by identifying important message flows. The authors use the concept of Shapley Value in Cooperative Game Theory and calculate Shapley Value as the initial assessments of flow importance score. A learning-based algorithm enables the refinement of important scores. Experiments are evaluated on both synthetic and real-world datasets.

With explainability in Graph Neural Networks (GNN) still in a nascent stage, most graph explanation methods generate explanations in terms of nodes, node features, edges, or sub-graphs. GNN's are essentially message-passing networks, where every node has access to a local view of the graph created by the propagation of neural messages (embeddings) along edges in the node’s local neighborhood. However, none of the previously proposed GNN explainability methods consider message-flows as possible explanations. The paper proposes FlowX to identify important message flows by employing concepts of Shapley values, where an approximation scheme is used to estimate the Shapley values as initial assessments of flow importance. Finally, a learning algorithm is proposed to refine and map the message flow importance to edge mask explanations, where experimental studies on both synthetic and real-world datasets show the improvement of explanations using FlowX.