L2XGNN: Learning to Explain Graph Neural Networks

Graph Neural Networks (GNNs) are a popular class of machine learning models. Inspired by the learning to explain (L2X) paradigm, we propose L2XGNN, a framework for explainable GNNs which provides faithful explanations by design. L2XGNN learns a mechanism for selecting explanatory subgraphs (motifs) which are exclusively used in the GNNs message-passing operations. L2XGNN is able to select, for each input graph, a subgraph with specific properties such as being sparse and connected. Imposing such constraints on the motifs often leads to more interpretable and effective explanations. Experiments on several datasets suggest that L2XGNN achieves the same classification accuracy as baseline methods using the entire input graph while ensuring that only the provided explanations are used to make predictions. Moreover, we show that L2XGNN is able to identify motifs responsible for the graph's properties it is intended to predict.


Introduction
Graph Neural Networks (GNNs) are a widely used class of machine learning models. Since graphs occur naturally in several domains such as chemistry, biology, and medicine, GNNs have experienced widespread adoption. Following a trend toward building more interpretable machine learning models, there have been numerous recent proposals to provide explanations for GNNs. Most of the existing approaches provide post-hoc explanations starting from an already trained GNN to identify edges and node attributes that explain the model's prediction. However, as highlighted in Faber et al. [2021], there might be some discrepancy between the ground-truth explanations and those attributed to the trained GNNs. Indeed, post-hoc explanations are often not able to faithfully represent the mechanisms of the original model Rudin [2018]. Unfortunately, the very definition of what constitutes a faithful explanation is still open to debate and there exist several competing positions on the matter. Recent work has also shown that post-hoc attribution methods are often not better than random baselines on the standard evaluation metrics for explanation accuracy and faithfulness Agarwal et al. [2022b]. Much fewer approaches have considered the problem of GNN explainability from an intrinsic perspective. In contrast to post-hoc methods, approaches with built-in interpretability provide explanations during training by introducing new mechanisms, e.g. prototypes , stochastic attention , or graph kernels Feng et al. [2022a]. Nonetheless, the introduction of new mechanisms to compute graph representations differ from standard GNNs computations. Therefore, the reasoning process of the above interpretable networks differ from the original GNN architectures making them not faithful by design. Our intent, instead, is to generate explanations for standard GNNs by keeping the computations as faithful as possible compared to the original network.
A recently proposed alternative to post-hoc methods is the learning to explain (L2X) paradigm . The core difference to post-hoc methods is that the models are trained to, in the forward pass, discretely select a small subset of the input features as well as the parameters of a downstream model that uses only the selected features to make a prediction. The selected features are, therefore, faithful by design as they are the only ones used by the downstream model. Since the subset of features is sampled discretely, L2X requires a method for computing gradients of an expectation over a discrete probability distribution.  proposed a gradient estimator based on a relaxation of the discrete samples and tailored to the k-subset distribution. However, since the original work only considers the case of selecting exactly-k features, directly applying prior methods to the graph learning tasks is not possible and requires significant changes.
With this work, we bring the L2X paradigm to graph representation learning. The important ingredient is a recently proposed method for computing gradients of an expectation over a complex exponential family distribution Niepert et al. [2021]. The method facilitates approximate gradient backpropagation for models combining continuously differentiable GNNs with a black-box solver of combinatorial problems defined on graphs. Crucially, this allows us to learn to sample subgraphs with beneficial properties such as being connected and sparse. Contrary to prior work, this also creates a dependency between the random variables representing the presence of edges. The proposed framework L2XGNN, therefore, learns to select explanatory subgraph motifs and uses these and only these motifs for its message-passing operations. To the best of our knowledge, this is the first method for learning to explain on standard GNNs. The proposed framework is extensible as it can work with any optimization algorithm for graphs imposing properties on the sampled subgraphs.
We compare two different sampling strategies for obtaining sparse subgraph explanations resulting from two optimization problems on graphs: (1) the maximum-weight k-edge subgraph and (2) the maximum-weight k-edge connected subgraph problem. We show empirically that L2XGNN when combined with a base GNN does not lose accuracy on several benchmark datasets. Moreover, we evaluate the explanations quantitatively and qualitatively. We also analyze the ability of L2XGNN to help in detecting shortcut learning which can be used for debugging the GNN. Given the characteristics of the proposed method, our work improves model interpretability and increases the clarity of known black-box models as GNNs while maintaining competitive predictive capabilities.

Background
Let G(V, E) be a graph with n = |V | the number of nodes. Let X ∈ R n×d be the feature matrix that associates each node of the graph with a d-dimensional feature vector and let A ∈ R n×n be the adjacency matrix. GNNs have three computations based on the message passing paradigm Hamilton et al. [2017] which is defined as where γ, , and φ represent update, aggregation and message function respectively. Propagation step. The message-passing network computes a message m ij = φ(h −1 i , h −1 j , r ij ) between every pair of nodes (v i , v j ). The function takes in input v i 's and v j 's representations h −1 i and h −1 j at the previous layer − 1, and the relation r ij between the two nodes. Aggregation step. For each node in the graph, the network performs an aggregation computation over the messages from The final embedding for node v i after L layers is z i = h L i and is used for node classification tasks. For graph classification, an additional readout function aggregates the node representations to obtain a graph representation h G . This function can be any permutation invariant function or a graph-level pooling function Ying et al. [2018], Zhang et al. [2018], Lee et al. [2019]. For Graph Isomorphism Networks (GINs) Xu et al. [2018], for instance, the message passing operation for node v i is where γ represents a multi-layer perceptron (MLP), and denotes a learnable parameter. We will write H = GNN (A, H −1 ) as a shorthand for the application of the th layer of the GNN under consideration.

Related Work
There are several methods to explain the behavior of GNNs. These methods use decomposition rules to decompose the model predictions leading back to the input space. The prediction is considered as the target score. Then, starting from the output layer, the target score is decomposed at each preceding layer according to the decided decomposition rules. In this way, the initial target score is distributed among the neurons at every layer. Finally, the decomposed terms obtained at the input layer are associated to the input features and used as importance scores of the corresponding nodes and edges.
Model-level methods. Yuan et al. [2020a]. Different from the instance-level methods above, these methods provide a general and high-level understanding of the models. In the context of GNNs, they aim at studying the input patterns that would lead to a certain target prediction. The generated explanations are general and provide a global understanding of the trained GNNs.
Prototype-based methods.  propose ProtGNN, a new explanatory method based on prototypes to provide built-in explanations, overcoming the limitations of post-hoc techniques. The explanations are obtained following case-based reasoning, where new instances are compared with several learned prototypes.
Concept-based methods. Magister et al. [2021] propose CGExplainer, a post-hoc explanatory methods for human-in-the-loop concept discovery. This concept representation learning method extracts concept-based explanations that allow the end-user to analyze predictions with a global view.
Additional works face the explainability problem from different perspectives as explanation supervision Gao et al. [2021], neuron analysis Xuanyuan et al. [2022], and motif-based generation Yu and Gao [2022]. For a comprehensive discussion on methods to explain GNNs, we refer the reader to the survey Yuan et al. [2020b]. A more detailed comparison with inherent interpretable methods and graph structure learning approaches can be found in the supplementary material.

Limitations of Prior Work
Existing XAI methods for GNNs have several limitations and can lead to inconsistencies Duval and Malliaros [2021], Faber et al. [2021]. We focus on the problem of identifying a subset of the edges as an explanation of the model's message-passing behavior. Hence, an explanation is equivalent to identifying a mask for the adjacency matrix of the original graph. Intuitively, an explanation can be accurate and/or faithful. It is accurate if it succeeds in identifying the edges in the input graph responsible for the graph's class label. This property can, for example, be evaluated with synthetic data where the class label of a graph is determined by the presence or absence of a particular substructure. An explanation is faithful if the edges identified as the explanation cause the prediction of the GNN on an input graph. Contrary to measuring accuracy, there is no consensus on evaluating faithfulness.
Recent work has proposed to measure unfaithfulness as the difference between the predictions of (1) the GNN on a perturbed adjacency matrix and (2)  , Agarwal et al. [2022b,a]. We believe that this definition is problematic as the perturbation is typically implemented using a swap operation which replaces two existing edges (a, b) and (c, d) with two new edges (a, c) and (b, d).
Hence, these new edges are present in the unmasked adjacency matrix but not present in the masked one. It is, however, natural that the same GNN would predict highly different label distributions on these two graphs. For instance, consider a chemical compound where we remove and add new bonds.
The resulting compounds and their properties can be chemically very different. Hence, contrary to prior work, we define a subgraph to be a faithful explanation, if it is a significantly smaller subgraph of the input graph and we know that only its structure is used in the message-passing operations of equation (1).

Learning to Explain Graph Neural Networks
We propose a method that learns both (i) the parameters of a graph generative model and (ii) the parameters of a GNN operating on sparse subgraphs approximately sampled from said generative model in the forward pass. In line with prior work on learning to explain , the maximum probability subgraph is then used at test time to make the prediction and, therefore, serves as the faithful explanation. Since we aim to sample graphs with certain properties (e.g., connected subgraphs) we need a new approach to sampling and gradient estimation. Contrary to prior work on edge masking Schlichtkrull et al. [2021] which treats edges as independent binary random variables, we use a recently introduced method for backpropagating through optimization algorithms. This allows us to select subgraphs with specific properties and, therefore, to explicitly model dependencies between edge variables.

Problem Statement and Framework
We aim to jointly learn the parameters of a probability distribution over subgraphs with certain properties and the parameters of a GNN operating on graphs sampled from said distribution in the context of the graph classification problem. Here, the training data consists of a set of triples {(A, X, y) j }, j ∈ {1, ..., N }, where A is an n × n binary adjacency matrix, X ∈ R n×d a node attribute matrix with d the number of node attributes, and y the target graph label. First, we have a learnable function h v : A × X → Θ where A is the set of all n × n adjacency matrices, X the set of all attribute matrices, v are the parameters of h, and Θ the set of possible edge parameter values. The function, which we refer to as the upstream model, maps the adjacency and attribute matrix to a matrix of edge weights θ ∈ R n×n . Intuitively, θ i,j is the prior probability of edge (i, j).
Next, we assume an algorithm opt : Θ → A which returns the (approximate) solutions to an optimization problem on edge-weighted graphs. Examples of such optimization problems are the maximum-weight spanning tree or the maximum-weight k-edge connected subgraph problems. The optimization algorithm is treated as a black box. One can choose the optimization problem according to the application's requirements. We have found, for instance, that the connected subgraphs lead to better explanations in the domain of chemical compound classification. Contrary to prior work, the optimization problem creates a dependency between the binary variables modeling the edges.
For every binary adjacency matrix Z ∈ A, we write Z ∈ F if and only if the adjacency matrix is a feasible solution (not necessarily an optimal one) of the chosen optimization problem. We can now define a discrete exponential family distribution as where ·, · F is the Frobenius inner product and B(θ) is the log-partition function defined as Hence, p is a probability distribution over adjacency matrices that are feasible solutions to the optimization problem under consideration. Each feasible adjacency matrix's probability mass is proportional to the product of its edge weights. For example, if the optimization problem is the maximum-weight k-edge connected subgraph problem, the distribution assigns a non-zero probability mass to all adjacency matrices of graphs that have k edges and are connected.
Given an optimization problem, we would like to sample exactly from the above probability distribution p(Z; θ). Unfortunately, this is intractable since computing the log-partition function is in general NP-hard. However, as in prior work Niepert et al. [2021], we can use perturb-and- MAP Papandreou and Yuille [2011] to approximately sample from the above distribution as follows. Let ∼ ρ( ) be a n × n matrix of appropriate random variables such as those following the Gumbel distribution. We can then approximately sample an adjacency matrix Z from p(Z; θ) by computing Z = opt(θ + ). Hence, we can approximately sample by perturbing the edge weights (unnormalized probabilities) θ and by applying the optimization algorithm to these perturbed weights.
In the final part of the model (the downstream model) we use the sampled Z as the input adjacency matrix to a message-passing neural network f u : In summary, we have the following model architecture for training input data (A, X, y): Figure 1 illustrates the architecture. With ω = (u, v) the learnable parameters of the model and the target variable y the loss is now defined as: with ] which can be estimated by Monte-Carlo sampling. In contrast, the gradient of L with respect to v is: is not continuously differentiable wrt θ. While it would be possible to use the score function estimator, its high variance makes it less competitive in practice Niepert et al. [2021].

Implicit Maximum-Likelihood Learning
The variant of I-MLE we use in this work estimates ∇ θ L(A, X, y; ω) by implicitly creating a target distribution q(Z; θ ) using perturbation-based implicit differentiation Domke [2010]. Here, the parameters θ are moved in the direction of −∇ Z (f u (A, X), y)), the negative gradient of the downstream loss with respect to the sampled adjacency matrix Z, to construct θ q(Z; θ ) := p(Z; θ − λ∇ Z (f u (Z, X), y)) (8) with Z = opt(θ + ) and λ > 0 the strength of the perturbation. Intuitively, by moving the weights θ into the direction of the negative gradients of Z, the resulting distribution q is more likely to generate samples with a lower downstream loss. We approximate ∇ θ L(A, X, y; ω) with Monte Carlo estimates of the gradients of the KL divergence between p and q: In other words, ∇ θ L(A, X, y; ω) is approximated by the difference between an approximate sample from p(Z; θ) and an approximate sample from q(Z; θ ). In this way we move the distribution p(Z; θ) closer to q(Z; θ ).
Algorithm 1 Greedy algorithm opt con for the maximum-weight k-edge connected subgraph problem.

L2XGNN: Learning to Explain GNNs with I-MLE
We now describe the class of L2XGNN models we use in the experiments. First, we need to define the function h v (A, X). Here we use a standard GNN (see equation 1) to compute for every node i and every layer the vector representation h i = h v (A, X) i,1:d . We then compute the matrix of edge weights by taking the inner product between each pair of node embeddings. More formally, we compute θ i,j = h i , h j for some fixed . Typically, we choose = 1.
In this work, we sample the noise perturbations from the sum of Gamma distribution Niepert et al. [2021]. Other noise distributions such as the Gumbel distribution are possible.

Sampling Constrained Subgraphs
An advantage of the proposed method is its ability to integrate any graph optimization problem as long as there exists an algorithm opt for computing (approximate) solutions. In this work, we focus on two optimization problems: (1) The maximum-weight k-edge subgraph and (2) the maximumweight k-edge connected subgraph problems. The former aims to find a maximum-weight subgraph with k edges. The latter aims to find a connected maximum-weight subgraph with k edges. Other optimization problems are possible but we found that sparse and connected subgraphs provide a good efficiency-effectiveness trade-off.
Computing maximum weight k-edge subgraphs is highly efficient as we only need to select the k edges with the maximum weights. In order to compute connected k-edge subgraphs we use a greedy approach. First, given a number k of edges, we select a single edge e i,j with the highest weight θ i,j from the input graph. At every iteration of the algorithm, we select the next edge such that it (a) is connected to a previously selected edge and (b) has the maximum weight among all those connected edges. A more detailed description of the greedy algorithm is given in Algorithm 1.
Finally, we need to define the function f u (the downstream function) of the proposed framework. Here, we again use a message-passing GNN that follows the update rule The neighborhood structure N (·), however, is defined through the output adjacency matrix Z of the optimization algorithm opt Hence, if after the subgraph sampling, there exists a node v i which is an isolated node in the adjacency matrix Z, that is, Z i,j = Z j,i = 0 ∀j ∈ {1, ..., n}, the embedding of the node will not be updated based on message passing steps with neighboring nodes. This means that, for isolated nodes, the only information used in the downstream model is the one from the nodes themselves. Conceptually, Z works as a mask over the messages m ij computed at each layer .
The adjacency matrix Z is then used in all subsequent layers of the GNN. In particular, for one layer we have where is the Hadamard product. Finally, the remaining part of the L2XGNN network for the graph classification is where we use a global pooling operator to generate the (sub)graph representation h G that will then be used by the MLP network η(·) to output a probability distributionŷ over the class labels. Finally, a loss function is applied whose gradients are used to perform backpropagation. At test time, we use the maximum-probability subgraph for the explanation and prediction, that is, we do not perturb at test time.

Experiments
First, we evaluate the predictive performance of the model compared to baselines. Second, we qualitatively and quantitatively analyze the explanatory subgraphs for datasets for which we know the ground-truth motifs. Finally, we analyze whether the generated output can be helpful for model debugging purposes. We report several ablation studies to investigate the effects of different model choices on the results in the supplementary material. The code for reproducing our experiments is available at https://github.com/GiuseppeSerra93/L2XGNN.

Datasets and Settings
Datasets. To understand the change in the predictive capabilities of the base models when integrating L2XGNN, we use six real-world datasets for graph classification tasks: MUTAG Debnath et al. [1991], PROTEINS Borgwardt et al. [2005], YEAST Yan et al. [2008], IMBD-BINARY, IMBD-MULTI Yanardag and Vishwanathan [2015], and DD Rossi and Ahmed [2015]. To quantitatively evaluate the quality of the explanations, we use datasets that include ground-truth edge masks. In particular, we use MUTAG 0 and BA2Motifs. MUTAG 0 is a dataset introduced in Tan et al. [2022] which contains the benzene-NO 2 (i.e., a carbon ring with a nitro group (NO 2 ) attached) as the only discriminative motif between positive and negative labels. A graphical representation of the benzene-NO 2 compound is given in Figure 2. BA2Motifs is a synthetic dataset that was first introduced in Luo et al. [2020]. The base graphs are Barabasi-Albert (BA) graphs. 50% of the graphs are augmented with a house-motif graphs, the rest with a 5-node cycle motif. The discriminative subgraph leading to different predictions is the motif attached to the BA graph.

Experimental settings.
To evaluate the quality of our approach, we use L2XGNN with several GNN base models including GCN Kipf and Welling [2017], GIN Xu et al. [2018] and GraphSAGE Hamilton et al. [2017]. We compare the results when using the original model and when the same model is combined with our XAI method. For model selection and evaluation, to fairly compare the methods, we follow a previously proposed protocol 1 . The hyperparameter selection is performed via 10-fold cross-validation, where a training fold is randomly used as a validation set. The selection is performed for the number of layers [1, 2, 3, 4] and the number of hidden units [16,32,64,128] based on the validation accuracy. For a fair comparison with the backbone architectures, we select the best configuration for each dataset, and we integrate our approach into the best model. Instead of fixing a value k for each input graph, we compute k based on a ratio R of edges to be kept. Once the hyperparameters of the default model are found, we select the best ratio R from the set of values [0.4, 0.5, 0.6, 0.7] based on the validation accuracy. We do not include extreme values for two reasons: (1) smaller values for R lead to reduced predictive capabilities and do not lead to meaningful explanatory subgraphs; and (2) higher values would not remove enough edges compared to the original input. Finally, we choose the perturbation intensity λ from the values [10, 100, 1000]. A summary of the selected hyperparameter values for each configuration can be found in the appendix.

Empirical Results
Graph Classification Comparison with Base GNNs. Following the experimental procedure proposed in , Table 1 lists the results of using L2XGNN with base GNN architectures  for graph classification tasks. The primary goal of this work is not to provide a better predictive model, but to provide faithful explanation masks while maintaining similar predictive performance. This analysis is important since inherent interpretable networks are known for creating a trade-off with the predictive capabilities of the model, and practitioners may not be willing to sacrifice the prediction accuracy for increased transparency . Nevertheless, we observe that L2XGNN is competitive and often even outperforms the base GNN models on the benchmark datasets. We train a 3-layer GIN for 200 epochs with hidden dimensions equal to 64 and a learning rate equal to 0.001. We save the best model according to the validation accuracy and we compare it with the post-hoc techniques. In our case, we integrate L2XGNN into the same architecture and learn the edge masking during training as described before. We report the graph classification results for the two datasets in the appendix. In Table 2, we report the explanation accuracy evaluation with respect to the ground-truth motifs in comparison with post-hoc techniques for 5 different data splits. The explanation problem is formalized as a binary classification problem, where the edges belonging to the ground-truth motif are treated as positive labels. We observe that L2XGNN obtains better or the same results as the considered explanatory models. While for the post-hoc explanation techniques we cannot guarantee that the GNNs use exclusively the explanation subgraphs for the prediction Yuan et al.
[2020b], our method, by providing faithful explanations, overcomes this limitation. It is exactly the provided explanation that is used in the message-passing operations of L2XGNN. Qualitative Evaluation of the Explanations. In Figure 3, we present some of the subgraphs identified by L2XGNN when combined with two different base GNNs. Based on prior studies and chemical domain knowledge Debnath et al. [1991], Lin et al. [2021], Tan et al. [2022], carbon rings (the black circles in the pictures) and NO 2 groups are known to be mutagenic. Interestingly, we can notice that, when using the information of connected subgraphs, the models are able to recognize a complete carbon ring with a NO 2 group in most of the cases. In some cases, the carbon ring is not complete, but the explanations are still helpful to understand which motifs are potentially important for the prediction. In the second row, we can observe the results of the sampling strategy where we GCN GIN w/ con w/o con Figure 3: Visualization of some of the subgraphs selected by L2XGNN for MUTAG0 on the test set. The solid edges are the ones sampled by our approach. The first and second rows depict, respectively, some motifs with or without constraining the subgraphs to be connected. Black, blue, red, and gray nodes represent carbon (C), nitrogen (N), oxygen (O), and hydrogen (H) atoms respectively.
do not require subgraphs to be connected. In this case, the carbon rings are not always identified. Instead, the NO 2 group is always considered important for the prediction. More generally, as also reported in Yuan et al. [2021], studying connected subgraphs results in more natural motifs compared to the motifs obtained without the connectedness constraint. The visualization of explanations for GraphSAGE and other methods can be found in the supplementary material.  Shortcut Learning Detection. By generating faithful subgraph explanations, our approach can be used to detect whether the predictive model is focusing on the expected features or if it is affected by shortcut learning. This is particularly important for GNNs, where seemingly small implementation differences can influence the learning process of the model Schlichtkrull et al. [2021]. To this end, we use the BA2Motifs dataset Luo et al. [2020]. We trained two different models, GCN and GIN, achieving a test accuracy of 0.67 and 1.0 respectively. Taking a closer look at the explanations of the first model, we observed that most of the correct predictions were (incorrectly) correlated with the cycle motif and that the explanations were similar to the ones reported in Figure 4. The explanatory results show that the model is not learning the expected discriminative motifs and, consequently, the accuracy for the test set is poor. This insight can help users to change the configuration of the architecture or to use a different model (e.g., GIN). More generally, the results highlight that faithful explanations can facilitate model analysis and debugging.

Conclusion and Limitations
We propose L2XGNN, a framework that can be integrated into GNN architectures to learn to generate explanatory subgraphs which are exclusively used for the models' predictions. Our experimental findings demonstrate that the integration of L2XGNN with base GNNs does not affect the predictive capabilities of the model for graph classification tasks. Furthermore, according to the definition provided in the paper, the resulting explanations are faithful since the retained information is the only one used by the model for prediction. Hence, differently from most of the common techniques, our explanations reveal the rationale of the GNNs and can also be used for model analysis and debugging.
A limitation of the approach is the reduced efficiency compared to baseline GNN models. Since we need to integrate an algorithm to compute (approximate) solutions to a combinatorial optimization problem, each forward-pass requires more time and resources. Moreover, depending on the choice of the optimization problem, we might not capture the structure of explanatory motifs required for the application under consideration.

Dataset Statistics
In Table 1, we report the statistics of the datasets used for graph classification tasks. We use datasets from different domains including biological data, social networks, and bioinformatics data. For a comprehensive evaluation, we include datasets with different characteristics, such as a larger number of graphs or a larger number of nodes and edges.

Comparison with Non-post-hoc Methods
Among the plethora of post-hoc methods for graphs, ProtGNN Zhang et al. [2021] and KerGNN Feng et al. [2022] are noteworthy exceptions. The first approach proposes a framework to generate explanations by comparing input graphs with prototypes learned during training, The latter, instead, combines graph kernels with the message passing paradigm to learn hidden graph filters. In contrast to these methods, our approach only relies on standard GNNs, without the use of prototypes or graph kernels. Consequently, even though both ProtGNN and KerGNN provide built-in explanations, they are not faithful by design as they introduce a new mechanism to compute graph representations which differs from standard GNNs computations. In terms of explanatory capability, the learned prototypes are not directly interpretable and need to be matched to the closest training subgraphs to be humanunderstandable. Graph filters, instead, do not necessarily match existing patterns in the instance-based case. In both cases, the output can only provide a general idea of the important structures used by the model for prediction but fail at revealing precisely the instance-level explanation for each input graph. For these reasons, a direct comparison with the two methods is not possible.

Comparison with Graph Structure Learning Approaches
Recently, there have been related methods for learning the structure of graph neural networks. Following the taxonomy proposed in Zhu et al. [2021], the structure learning methods most related to L2XGNN fall into the postprocessing category, and more specifically, under the discrete sampling subcategory. All existing methods use variants of the Gumbel-softmax trick which is limited in modeling complex distributions. Moreover, only when the straight-through version of the Gumbelsoftmax trick is used, one can obtain truly discrete and not merely relaxed adjacency matrices in the forward pass. In contrast, L2XGNN always samples purely discrete adjacency matrices. It is, to the best of our knowledge, the only method that allows us to model complex dependencies between the edge variables through its ability to integrate a combinatorial optimization algorithm on graphs. Other strategies include sampling edges between each pair of nodes from a Bernoulli distribution Franceschi et al. [2019] or sampling subgraphs for subgraph aggregation methods in a data-driven manner Qian et al. [2022]. All these methods, however, are not concerned with the problem of explaining the behavior of GNNs explicitly.

Ablation Study
In the graph classification task, we compare the two sampling strategies. From the results, the connected sampling is able to get better results than the non-connected counterpart on most datasets. In fact, the connectivity of subgraphs is essential to grasp the complete information about the important patterns, especially for chemical compound data where connected atoms are usually expected to create molecules or chemical groups. This aspect is also supported by the results obtained in the explanation accuracy task, where the connected strategy returns better explanations for the chemical dataset. Additionally, as previously mentioned, evaluating connected structures rather than just important edges looks more natural and intelligible. In Table 3, we analyze the effect of the quantity of retained information on the prediction accuracy. A smaller ratio indicates that we retain fewer edges during training and, consequently, the resulting subgraphs are more sparse and, therefore, interpretable. As one can see, this affects the predictive capabilities only when R is small. Starting from R = 0.5, the ratio does not affect particularly the predictive capabilities of the model. In fact, for graph classification tasks, some of the information contained in the initial computational graph does not condition the prediction as the information may be redundant or noisy. For instance, considering the MUTAG dataset, we know that the initial graphs contain on average 20 edges. The discriminative motif benzene-NO 2 , instead, contains around 9 edges, meaning that we ideally need 50% of the original edges to obtain good results. This is in line with the findings of this analysis and the graph classification results previously reported in the main paper.

Time Complexity Analysis
As reported in Feng et al. [2022], most message-passing architectures have a time complexity of O(n 2 ). Thus, since L2XGNN uses native GNNs, it also has a worst-case complexity of O(n 2 ). KerGNN instead, being based on graph kernels, has a time complexity between O(n 2 ) and O(n 3 ) in case the graph is fully-connected. The time overhead of opt depends on the algorithm used. For the maximum-weight k-edge subgraph problem, we only need to select the k edges with the maximum weights. To find the maximum-weight k-edge connected subgraph, the greedy algorithm needs to find the maximum weight edge among all edges adjacent to the currently selected subgraph which, in the worst case, are all edges. Finally, to conclude the comparison with post-hoc techniques, we analyze the time efficiency for generating explanations for unseen data. For a coherent comparison, we use the same test splits used for the previous experiments. Table 4 reports the average time cost to obtain the explanations related to the graphs in the test sets.

Graph Classification Accuracy for Synthetic Datasets
In Table 5, we report the graph classification accuracy for the datasets used in our explanation accuracy experiment. In particular, we compare the 3-layer GIN architecture used for generating explanations through post-hoc techniques and the same architecture when our approach is integrated during training. The data splits are the same used in the previous evaluation. Again, the results demonstrate that the integration of L2XGNN does not degrade significantly the predictive capabilities of the original model.

Visual Comparison of Generated Explanations
In Figure 1 we provide a visual analysis of the explanations generated by the methods considered in our experiments. The graph visualization supports the numerical evaluation carried on in the main paper. In fact, one can see that the explanations generated by the post-hoc approaches may vary substantially depending on the given input graph. In our case, instead, the explanations remain constant regardless of the input information. This claim supports the explanation accuracy analysis, where our approach has one of the smallest standard deviation among all the considered methods. Additionally, we also included the explanations generated with an attention-based GNN, namely GAT Veličković et al. [2018]. Although having a similar predictive performance in the graph classification task (99.6 ± 0.03), the resulting explanations are not qualitatively comparable with our approach. This is in line with previous work , ,  asserting that graph attention models are not able to generate attention weights with high-fidelity, and consequently, cannot provide faithful and meaningful explanations.

GAT
GraphSAGE dsc GraphSAGE con Figure 1: Comparison of the generated explanations for MUTAG0 on the test set. The solid edges are the ones considered responsible of a correct prediction. Black, blue, red, gray, and green nodes represent carbon (C), nitrogen (N), oxygen (O), hydrogen (H), and chlorine (Cl) atoms respectively.