Graph Classification via Reference Distribution Learning: Theory and Practice

Authors: zixiao  wang
jicong  fan
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 15:08:51

Abstract: Graph classification is a challenging problem owing to the difficulty in quantifying the similarity between graphs or representing graphs as vectors, though there have been a few methods using graph kernels or graph neural networks (GNNs). Graph kernels often suffer from computational costs and manual feature engineering, while GNNs commonly utilize global pooling operations, risking the loss of structural or semantic information. This work introduces Graph Reference Distribution Learning (GRDL), an efficient and accurate graph classification method. GRDL treats each graph's latent node embeddings given by GNN layers as a discrete distribution, enabling direct classification without global pooling, based on maximum mean discrepancy to adaptively learned reference distributions. To fully understand this new model (the existing theories do not apply) and guide its configuration (e.g., network architecture, references' sizes, number, and regularization) for practical use, we derive generalization error bounds for GRDL and verify them numerically. More importantly, our theoretical and numerical results both show that GRDL has a stronger generalization ability than GNNs with global pooling operations. Experiments on moderate-scale and large-scale graph datasets show the superiority of GRDL over the state-of-the-art, emphasizing its remarkable efficiency, being at least 10 times faster than leading competitors in both training and inference stages.

Main paper: 
 introduction: 
 Graphs serve as versatile models across diverse domains, such as social networks [Wang et al., 2018] , biological compounds [Jumper et al., 2021] , and the brain [Ktena et al., 2017] . There has been considerable interest in developing learning algorithms for graphs, such as graph kernels [Gärtner et al., 2003 , Shervashidze et al., 2011 , Chen et al., 2022b] and graph neural networks (GNNs) [Kipf and Welling, 2016 , Defferrard et al., 2016 , Gilmer et al., 2017] . GNNs have emerged as powerful tools, showcasing state-of-the-art performance in various graph prediction tasks [Veličković et al., 2017 , Gilmer et al., 2017 , Hamilton et al., 2017 , Xu et al., 2018 , Sun et al., 2019 , You et al., 2021 , Ying et al., 2021 , Liu et al., 2022b , Chen et al., 2022a , Xiao et al., 2022 , Sun et al., 2023 , Sun and Fan, 2024 , Sun et al., 2024] . Despite the evident success of GNNs in numerous graph-related applications, their potential remains underutilized, particularly in the domain of graph-level classification.Current GNNs designed for graph classification commonly consist of two components: the embedding of node features through message passing [Gilmer et al., 2017] and subsequent aggregation by some permutation invariant global pooling (also called readout) operations [Xu et al., 2018] . The primary purpose of pooling is to transform a graph's node embeddings, a matrix, into a single vector.Empirically, pooling operations play a crucial role in classification [Ying et al., 2018] . However, these pooling operations tend to be naive, often employing methods such as simple summation or averaging. These functions collect only first-order statistics, leading to a loss of structural or semantic information. In addition to the conventional sum or average pooling, more sophisticated pooling operations have shown improvements in graph classification [Li et al., 2015 , Ying et al., 2018 , Lee et al., 2019 , 2021 , Buterez et al., 2022 , Yu et al., 2024] , but they still carry the inherent risk of information loss.Different from graph kernel methods and existing GNN methods, we propose a novel GNN method that classifies the nodes' embeddings themselves directly, thus avoiding the global pooling step. In our method, we treat the nodes' latent representations of each graph, learned by a neural network, as a discrete distribution and classify these distributions into K different classes. The classification is conducted via measuring the similarity between the latent graph's distributions and K discriminative reference discrete distributions. The reference distributions can be understood as nodes' embeddings of representative virtual graphs from K different classes, and they are jointly learned with the parameters of the neural network in an end-to-end manner. To evaluate our method, we analyze the generalization ability of our model both theoretically and empirically. Our contributions are two-fold.• We propose a novel graph classification method GRDL that is efficient and accurate.-GRDL does not require any global pooling operation and hence effectively preserves the information of node embeddings. -Besides its high classification accuracy, GRDL is scalable to large graph datasets and is at least ten times faster than leading competitors in both training and inference stages. • We provide theoretical guarantees, e.g. generalization error bounds, for GRDL.-The result offers valuable insights into how the model performance scales with the properties of graphs, neural network structure, and reference distributions, guiding the model design. -For instance, the generalization bounds reveal that the references' norms and numbers have tiny impacts on the generalization, which is also verified by the experiments. -More importantly, we theoretically prove that GRDL has a stronger generalization ability than GNNs with global pooling operations.The rest of this paper is organized as follows. We introduce our model in Section 2 and analyze the generalization ability in Section 3. Related works are discussed in Section 4. Section 5 presents the numerical results on 11 benchmarks in comparison to 12 competitors.2 Proposed Approach 
model framework: 
 Following convention, we denote a graph with index i by G i = (V i , E i ), where V i and E i are the vertex (node) set and edge set respectively. Given a graph dataset G = {(G 1 , y 1 ), (G 2 , y 2 ), . . . , (G N , y N )}, where y i ∈ {1, 2, . . . , K} is the associated label of G i and y i = k means G i belongs to class k, the goal is to learn a classifier f from G that generalizes well to unseen graphs. Since in many scenarios, each node of a graph has a feature vector x and the graph is often represented by an adjacency matrix A, we also write G i = (A i , X i ) for convenience, whereA i ∈ R ni×ni , X i ∈ R ni×d0 , n i = |V i |is the number of nodes of graph i, and d 0 denotes the number of features. We may alternatively denote the graph dataset as G = {((A 1 , X 1 ), y 1 ), ((A 2 , X 2 ), y 2 ), . . . , ((A N , X N ), y N )}.Our approach is illustrated in Figure 1 . For graph classification, we first use a GNN, denoted as f G , to transform each graph to a node embedding matrix H i ∈ R ni×d that encodes its properties, i.e.,EQUATIONwhere f G ∈ F G and F G denotes a hypothesis space. The remaining task is to classify H i without global pooling. Direct classification of node embeddings is difficult due to two reasons:(i) Different graphs have different numbers of nodes, i.e. in general, n i ̸ = n j if i ̸ = j.(ii) The node embeddings of each graph are permutation invariant, namely, PH i and H i represent the same graph for any permutation matrix P.However, the two properties are naturally satisfied if we treat the node embeddings of each graph as a discrete distribution. Specifically, each H i is a discrete distribution and each row of H i is anFigure 1 : The GRDL framework. Classification involves using a GNN f G to encode a graph's information into a node embedding distribution. The similarities between the node embeddings and K reference distributions are calculated by the reference module f D . The graph is assigned the label of the reference that exhibits the highest similarity.outcome of the distribution. There is no order between the outcomes in each distribution. Also, different distributions may have different numbers of outcomes. Before introducing our method in detail, we first give a toy example where the commonly used mean and max pooling operations fail.Example 2.1. Suppose two graphs G 1 and G 2 have self-looped adjacency matrices Ã1 = [1 1 1 0; 1 1 0 1; 1 0 1 1; 0 1 1 1] and Ã2 = [1 0 0 0; 0 1 1 0; 0 1 1 0; 0 0 0 1] respectively, and have one-dimensional node features X 1 = [3 6 9 12] ⊤ and X 2 = [6 6 9 9] ⊤ respectively. Let Âi be the normalized adjacency matrices, i.e., Âi = diag( Ãi 1) -1/2 Ãi diag( Ãi 1) -1/2 . Performing the neighbor aggregation H i = Âi X i , i = 1, 2, we obtain H 1 = [6 7 8 9] ⊤ and H 2 = [6 7.5 7.5 9] ⊤ . We see that mean(H 1 ) = mean(H 2 ) = 7.5 and max(H 1 ) = max(H 2 ) = 9. This means the simple mean and max pooling operations failed to distinguish the two graphs. In contrast, our method treats H 1 and H 2 as two different discrete distributions and hence is able to distinguish the two graphs. Note that incorporating a learnable parameter W, i.e., H i = Âi X i W, or performing multiple times of neighbor aggregation does not change the conclusion.We propose to classify the discrete distributions {H 1 , H 2 , . . . , H N } ≜ H by a reference layer f D .The classification involves measuring the similarity between H i and K reference discrete distributions {D 1 , D 2 , . . . , D K } ≜ D that are discriminative. Each D k ∈ R m k ×d can be understood as node embeddings of a virtual graph from the k-th class, k ∈[K]. We make m 1 = • • • = m K = m for convenience.Letting ξ be a similarity measure between two discrete distributions, thenEQUATIONThis forms a matrix S = [s 1 , s 2 , . . . , s N ] ⊤ whereEQUATIONf D ∈ F D and F D denotes a hypothesis space induced by the reference layer. References in D are parameters of the reference layer and are jointly learned with node-embedding network parameters in an end-to-end manner. Now combining Equation (1) and Equation (3), we arrive atEQUATIONf D • f G calculates s i , representing similarities between G i and all references. We get G i 's label byEQUATIONTo train the model, we first use the softmax function to convert s i to a label vector ŷi = [ŷ i1 , . . . , ŷiK ] ⊤ , whereEQUATIONUsing the cross-entropy loss, we minimizeL CE = - 1 N N i=1 K k=1 y ik log ŷik . (7)Intuitively, the reference distributions in D should be different from each other to ensure discriminativeness. Therefore, we also consider the following discrimination loss:EQUATIONThen we solve the following problem:EQUATIONwhere λ ≥ 0 is a hyperparameter. We call (9) Graph Classification via Reference Distribution Learning (GRDL). Specific designs of F G and F D are detained in the next section.We get GRDL's network F by concatenating the node embedding module and the reference module:EQUATION)Design of F G We use an L-layer message passing network as our node embedding module F G :EQUATION)F l is the l-th message passing layer (e.g. a GIN layer [Xu et al., 2018] ) that updates the representation of a node by aggregating representations of its neighbors, meaningEQUATION)where h(l) v is the feature vector of node v produced by the l-th layer F l . Different GNNs have different choices of COMBINE (l) (•) and AGGREGATE (l) (•).Design of F D Based on (3), the hypothesis space defined by the reference layer isEQUATIONIn our work, we choose ξ(•, •) to be the negative squared Maximum Mean Discrepancy (MMD).Initially used for two-sample tests, MMD is now widely used to measure the dissimilarity between distributions [Gretton et al., 2012a] . For an embedding H ∈ R n×d and a referenceEQUATION)where ϕ is some feature map, h ⊤ i is the i-th row of H, and d ⊤ j is the j-th row of D. The MMD in ( 14) is known as biased MMD [Gretton et al., 2012a] and its performance is almost the same as the unbiased one, in our experiments. Therefore we only present (14) here. Using kernel trick k(x, x ′ ) = ϕ(x) ⊤ ϕ(x ′ ), we obtain from (14) that ξ(H, D) = 2 mn n i=1 m j=1 k(h i , d j ) - 1 n 2 n i=1 n i ′ =1 k(h i , h i ′ ) - 1 m 2 m j=1 m j ′ =1 k(d j , d j ′ ).In this work, we employ the Gaussian kernel, i.e.,k(x, x ′ ) = exp -θ∥x -x ′ ∥ 2 2 (15)where θ > 0 is a hyperparameter. The Gaussian kernel defines an infinite-order polynomial feature map ϕ, covering all orders of statistics of the input variable. Consequently, MMD with the Gaussian kernel characterizes the difference between two distributions across all moments. Actually, we found that, in GRDL, the Gaussian kernel often outperformed other kernels such as the polynomial kernel.Several other statistical distances are available for measuring the difference between distributions, including Wasserstein distance and Sinkhorn divergence [Peyré and Cuturi, 2020] . However, their computational complexity is prohibitively high, making the model impractical for large-scale graph datasets. We also find, through experiments, that in our method, the classification performance of MMD is better than that of Wasserstein distance and Sinkhorn divergence as shown later in Table 1 . These explain why we prefer MMD.The θ in the Gaussian kernel (15) plays a crucial role in determining the statistical efficiency of MMD.Optimally setting of θ remains an open problem and many heuristics are available [Gretton et al., 2012b] . To simplify the process, we make θ learnable in our GRDL and rewrite ξ as ξ θ . Our empirical results in Appendix D.5 show that GRDL with learnable θ performs better. For convenience, we denote all the parameters of f G as w and let f w,D,θ = f D • f G . Then we rewrite problem (9) asEQUATIONThe (mini-batch) training of GRDL model is detailed in Algorithm 1 (see Appendix C). 
theoretical analysis: 
 In this section, we provide theoretical guarantees for GRDL, due to the following motivations:• As the proposed approach is novel, it is necessary to understand it thoroughly using theoretical analysis, e.g., understand the influences of data and model properties on the classification. • It is also necessary to provide guidance for the model design to guarantee high accuracy in inference stages.Matrix constructions We construct big matrices X, A and D, whereX = X ⊤ 1 , X ⊤ 2 , . . . , X ⊤ N ⊤ ∈ R ( i ni)×d ; A = diag(A 1 , A 2 , . . . , A N ) ∈ R ( i ni)×( i ni) is a block diagonal matrix, D = D ⊤ 1 , D ⊤ 2 , . . . , D ⊤ K ⊤ ∈ R Km×d .The adjacency matrix with self-connectivity is Ã = A + I. The huge constructed graph is denoted by G = ( Ã, X). This construction allows us to treat all graphs in dataset G as a whole and it is crucial for our derivation.Neural network Previously, for a deterministic network f ∈ F, its output after feeding forward a single graph is f (G i ). However, we mainly deal with the huge constructed graph G in this section, and notation will be overloaded tof (G) = S ∈ R N ×K , a matrix whose i-th row is f (G i ) ⊤ .We instantiate the message passing network as Graph Isomorphism Network (GIN) [Xu et al., 2018] . We choose to focus on GIN for two reasons. Firstly, the analysis on GIN is currently limited, most of the current bounds for GNNs don't apply for GIN [Garg et al., 2020 , Liao et al., 2021 , Tang and Liu, 2023] . The other reason is that GIN is used as the message-passing network in our numerical experiments. Notably, our proof can be easily adapted to other message-passing GNNs (e.g. GCN [Kipf and Welling, 2016] ). GIN updates node representations asEQUATION)where h(l)v denotes the node features generated by l-th GIN message passing layer. Let ε (l) = 0 for all layers and suppose all MLPs have r layers, the node updates can be written in matrix form asEQUATION)whereW (l) i ∈ R d (l) i-1 ×d (l) iis the weight matrix, and H (l) is the matrix of node features with H (0) = X. σ(•) is the non-linear activation function. Let F l be the function space induced by the l-th message passing layer, meaningEQUATIONwhere B(l)i is some constraint set on the weight matrix W (l)i and H (l) is given by ( 18). The L-layer GIN function space F G is the composition of F l for l ∈ [L], i.e.,EQUATION)Lettings ik = -MMD 2 (H (L)i , D k ), the reference layer defines the following function spaceEQUATION)Our proposed network (GRDL) is essentiallyF := F D • F G .Loss Function Instead of the cross entropy loss (7), we consider a general loss function l γ (•, •) satisfying 0 ≤ l γ ≤ γ to quantify the model performance. Importantly, this loss function is not restricted to the training loss because our generalization bound is optimization-independent. For instance, the loss function can be the ramp loss that is commonly used for classification tasks [Bartlett et al., 2017 , Mohri et al., 2018] . Given a neural network f ∈ F, we want to upper bound the model population risk of graphs and labels from an unknown distribution X × YEQUATIONGiven the observed graph dataset G sampled from X × Y, the empirical risk isEQUATIONof which ( 7) is just a special case. Appendix E provides more details about the setup and our idea.For convenience, similar to [Bartlett et al., 2017 , Ju et al., 2023] , we make the following assumptions. Assumption 3.1. The following conditions hold for F γ := {(G, y) → l γ (f (G), y) : f ∈ F}:(i) The activation function σ(•) is 1-Lipschitz (e.g. Sigmoid, ReLU).(ii) The weight matrices satisfyW (l) i ∈ B (l) i := {W (l) i : ∥W (l) i ∥ σ ≤ κ (l) i , ∥W (l) i ∥ 2,1 ≤ b (l) i }. (iii) The constructed reference matrix satisfy ∥D∥ 2 ≤ b D . (iv) The Gaussian kernel parameter θ is fixed. (v) The loss function l γ (•, y) : R K → R is µ-Lipschitz w.r.t ∥ • ∥ 2 and 0 ≤ l γ ≤ γ. Theorem 3.2 (Generalization bound of GRDL). Let n = min i n i , c = ∥ Ã∥ σ , and d = max i,l d (l) i . Denote R G := c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 ( r i=1 κ (l) i ) 2 L l=1 r i=1 b (l) i κ (l) i 2/3 3 . For graphs G = {(G i , y i )} N i=1drawn i.i.d from any probability distribution over X × {1, . . . , K} and references{D k } K k=1 , D k ∈ R m×d, with probability at least 1 -δ, every loss function l γ and network f ∈ F under Assumption 3.1 satisfyL γ (f ) ≤ Lγ (f ) + 3γ ln (2/δ) 2N + 8γ+24 √ v1+v2 ln N +24γ √ N v2 ln v3 N where v 1 = 64θKR G µ 2 n , v 2 = Km d, and v 3 = 24 √ θN b D µ √ m.The bound shows how the properties of the neural network, graphs, reference distributions, etc, influence the gap between training error and testing error. A detailed discussion will be presented in Section 3.3. Some interesting corollaries of Theorem 3.2, e.g., misclassification rate bound, can be found in Appendix F.7. Besides small generalization error L γ (f ) -Lγ (f ), a good model should have small empirical risk Lγ (f ). The empirical risk Lγ (f ) is typically a surrogate loss of misclassification rate of training data and a lower misclassification rate implies a smaller Lγ (f ). We now provide a guarantee for the correct classification of training data, namely small Lγ (f ).Notably, the node embeddings H i from the k-th class as well as the reference distributions D k are essentially some finite samples from an underlying continuous distribution P k . One potential risk is that, although the continuous distributions P 1 , P 2 , . . . , P K are distinct, we can only observe their finite samples and may fail to distinguish them from each other with MMD. Specifically, suppose a node embedding H i is from the k-th class, although 0 = MMD(P k , P k ) < MMD(P k , P j ) for any j ̸ = k, it is likely that MMD(H i , D k ) > MMD(H i , D j ) for some j ̸ = k. The following theorem provides the correctness guarantee for the training dataset G:Theorem 3.3. All graphs in the training set G are classified correctly with probability at least 1 -δ ifmin i̸ =j MMD(P i , P j ) > 1 √ m + 1 √ n 4 + 4 log 2N δ .Theorem 3.3 implies that a larger reference distribution size m benefits the classification accuracy of training data, resulting in a lower Lγ (f ). Moreover, a larger min i̸ =j MMD(P i , P j ) also makes correct classification easier according to the theorem, justifying our usage of discriminative loss (8).Let κ = max i,l κ(l) i , b = max i,l b (l) i κ (l) iand suppose δ is large enough, we simplify Theorem 3.2 asL γ (f ) ≤ Lγ (f ) + Õ √ v1+γ √ N v2 N ≤ Lγ (f ) + Õ µ b∥X∥2c L (Lr) 3 2 κLr √ θK/n N + γ Kmd N I.Dependence on graph property One distinctive feature of our bound is its dependence on the spectral norm of graphs' adjacency matrix. The large adjacency matrix Ã is a block-diagonal matrix, so its spectral norm c = ∥ Ã∥ σ = max i∈[N ] ∥ Ãi ∥ σ . By Lemma F.8, incorporating c L is sufficient for any L-step GIN message passing. This result aligns with Ju et al. [2023] , who achieved this conclusion via PAC-Bayesian analysis. Our derivation, based on the Rademacher complexity, provides an alternative perspective supporting this result. Notably, Liao et al. [2021] and Garg et al. [2020] proposed bounds scaling with graphs' maximum node degree, which is larger than the spectral norm of the graphs' adjacency matrix (Lemma F.18). Consequently, our bound is tighter.II. Use moderate-size message passing GIN The bound scales with the size of the message passing GIN, following Õ(c L (Lr) 3 2 κLr ). Empirical observations reveal κ > 1, and we prove that c > 1 (refer to Lemma F.20). Therefore, when the message-passing GNN has sufficient expressive power (resulting in a small Lγ (f )), a network with a smaller L and r may guarantee a tighter bound on the population risk compared to a larger one. Therefore, a promising strategy is to use a moderate-size message passing GNN. This is empirically supported by Figure 5 of Appendix D.7.The bound scales with the size of reference distributions m as Õ( √ m). When m is smaller, the bound tends to be tighter. However, if m is too small, the model's expressive capacity is limited, potentially resulting in a large empirical risk Lγ (f ), and consequently, a large population risk. Therefore, using moderate-size references is a promising choice, as supported by our empirical validation results in Appendix D.3 (see Figure 6 ).IV. Regularization on references norm barely helps Regularizing the norm of references ∥D∥ 2 , i.e., reducing b D , might be considered to enhance the model's generalization. However, it is important to note that b D only influences the term v 3 (in logarithm) in Theorem 3.2 and has a tiny influence on the overall bound. Conversely, such regularization constrains the model's expressive capacity, potentially leading to a large Lγ (f ) and increasing the population risk. This observation is empirically supported by experiments in Appendix D.7 (see Table 10 ).V. GRDL has a tighter bound than GIN with global pooling In Appendix A, we provide the generalization error bound, i.e., Theorem A.1, for GIN with global pooling and compare it with Theorem 3.2. The result shows that our GRDL has a stronger generalization ability than GIN, which is further supported by the numerical results in Table 4 . Remark 3.4. Currently, we use K reference distributions for classification (one for each class). One natural approach to enhancing the model's expressive power is increasing the number of references for each class. However, counterintuitively, our empirical observations, supported by Theorem B.1, suggest that having only one reference per class is optimal. We discuss this further in Appendix B. 
related work: 
 Various sophisticated pooling operations have been designed to preserve the structural information of graphs [Bianchi et al., 2020 , Ranjan et al., 2020 , Baek et al., 2021 , Chen and Gel, 2023 , Yu et al., 2024] . For instance, DIFFPOOL, designed by Ying et al. [2018] , learns a differentiable soft cluster assignment for nodes and maps nodes to a set of clusters to output a coarsened graph. Another method by Lee et al. [2019] utilizes a self-attention mechanism to distinguish nodes for retention or removal, and both node features and graph topology are considered with the self-attention mechanism.A recent research direction focuses on preserving structural information by leveraging the optimal transport (OT) [Peyré and Cuturi, 2020] . OT-GNN, proposed by Chen et al. [2021] , embeds a graph to a vector by computing Wasserstein distances between node embeddings and some "learned point clouds". TFGW, introduced by Vincent-Cuaz et al. [2022] , embeds a graph to a vector of Fused Gromov-Wasserstein (FGW) distance [Vayer et al., 2018] to a set of "template graphs". OT distances have also been combined with dictionary learning to learn graph vector embedding in an unsupervised way (GDL) [Liu et al., 2022a , Vincent-Cuaz et al., 2021 , Zeng et al., 2023] .Similar to the "learned point clouds" in OT-GNN, "template graphs" in TFGW, and dictionaries in GDL, our GRDL preserves information in node embeddings using reference distributions. To the best of the authors' knowledge, we are the first to model a graph's node embeddings as a discrete distribution and propose to classify it directly without aggregating it into a vector, marking our novel contribution. Additionally, our work stands out as the first to analyze the generalization bounds for this type of model, adding a theoretically grounded dimension to the research. By the way, our method is much more efficient than OT-GNN and TFGW. Please see Figure 2 and Table 8 .5 Numerical ExperimentsQuestion: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)? Answer: [Yes] Justification: The paper fully discloses all the information needed to reproduce the main experimental results. Guidelines:• The answer NA means that the paper does not include experiments.• If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not. • If the contribution is a dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable. • Depending on the contribution, reproducibility can be accomplished in various ways.For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, it may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed. • While NeurIPS does not require releasing code, the conference does require all submissions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. (d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome to describe the particular way they provide for reproducibility.In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results. 
graph classification benchmark: 
 Datasets We leverage eight popular graph classification benchmarks [Morris et al., 2020] , comprising five bioinformatics datasets (MUTAG, PROTEINS, NCI1, PTC-MR, BZR) and three social network datasets (IMDB-B, IMDB-M, COLLAB). We also use three large-scale imbalanced datasets (PC-3, MCF-7, and ogbg-molhiv [Hu et al., 2020] ). A summary of data statistics is in Table 6 . Baselines Our approach is benchmarked against four groups of state-of-the-art baselines: 1) GNN models with global or sophisticated pooling operations, including PATCHY-SAN [Niepert et al., 2016] , DIFF-POOL [Ying et al., 2018] , GIN [Xu et al., 2018] , DropGIN [Papp et al., 2021] , SEP [Wu et al., 2022] , GMT [Baek et al., 2021] , MinCutPool [Bianchi et al., 2020] , ASAP [Ranjan et al., 2020] , and Wit-TopoPool [Chen and Gel, 2023] ; 2) Optimal transport based models such as WEGL [Kolouri et al., 2020] and OT-GNN [Chen et al., 2021] ; 3) Kernel-based approaches including FGW [Titouan et al., 2019] operating on adjacency (ADJ) and shortest path (SP) matrices, the WL subtree kernel [Shervashidze et al., 2011] , and the Wasserstein WL kernel [Togninalli et al., 2019] ; 4) Graph transformers including Graphormer [Ying et al., 2021] and SAT [Chen et al., 2022a] . We also show the results of two variations of our GRDL: GRDL using Sinkhorn divergence (GRDL-S) and GRDL using Wasserstein distance (GRDL-W). For large imbalanced datasets, we only benchmark our GRDL against PATCHY-SAN, GIN, and DIFFPOOL because other methods are too costly. Details about the initialization and hyper-parameters setting can be found in Appendix D.3.Experiment Settings Due to the page limitation, please refer to Appendix D.2.Table 1 shows the classification results. The AUC-ROC scores of experiments results on the three large imbalanced datasets are reported in Table 2 . Our method has top 3 classification performance over baselines in almost all datasets. Our GRDL, GRDL-W and GRDL-S have close performance. However, as shown later in Figure 2 , our original GRDL has significantly lower time costs and thus is preferable for practical use. Graph transformers also have competitive performance, but they have significantly larger amount parameters and much higher time costs than our model, as shown by Table 13 in Appendix D.8 .We compare the time cost of our GRDL with two models that leverage optimal transport distances discussed in Section 4: OT-GNN [Chen et al., 2021] and TFGW [Vincent-Cuaz et al., 2022] .Compared with them, our model has significantly lower time costs. We present empirical average training time per epoch in Figure 2 and average prediction time per graph in Table 9 in Appendix D.4. Experiments were conducted on CPUs (Apple M1) using identical batch sizes, ensuring a fair comparison. It's noteworthy that the OT solver employed in TFGW and OT-GNN is currently confined to CPU, influencing the choice of hardware for this evaluation. We analyzed the theoretical time complexity in Appendix D.4 (see Table 8 ). We also compare training time with two latest pooling methods including Wit-TopoPool [Chen and Gel, 2023] and MSGNN [Lv et al., 2023] on eight real datasets and three synthetic datasets. The three synthetic datasets have 2000 graphs with 100(SYN-100), 300(SYN-300), and 500(SYN-500) nodes per graph, respectively. The edge number is 0.1n 2 where n is the number of nodes. The empirical training time per epoch is shown in Table 3 , where empty of MSGNN means it takes more than 200 seconds to train a single epoch, which is too costly. As can be seen, our method is the most efficient among these three methods.We use t-SNE [Van der Maaten and Hinton, 2008] to visualize the distributions of graphs' node embeddings given by our GRDL model, which is equivalent to visualizing each graph in a 3-D coordinate system. Firstly we use MMD to calculate a distance matrix C ∈ R (N +K)×(N +K) between G R D L G R D L -W G R D L -S O T -G N N T F G W 0 2 4 6 8 Time Cost (second) 0.4 1.8 2.5 3.6 9.3 MUTAG G R D L G R D L -W G R D L -S O T -G N N T F G WQuestion: Does the paper provide open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results, as described in supplemental material? Answer: [Yes] Justification:The paper provides open access to the data and code, with sufficient instructions to faithfully reproduce the main experimental results. Guidelines:• The answer NA means that paper does not include experiments requiring code.• Please see the NeurIPS code and data submission guidelines (https://nips.cc/ public/guides/CodeSubmissionPolicy) for more details. • While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
conclusions: 
 We proposed GRDL, a novel framework for graph classification without global pooling operations and hence effectively preserve the information of node embeddings. What's more, we theoretically analyzed the generalization ability of GRDL, which provided valuable insights into how the generalization ability scales with the properties of the graph data and network structure. Extensive experiments on moderatescale and large-scale benchmark datasets verify the effectiveness and efficiency of GRDL in comparison to baselines. However, on some benchmark datasets (e.g. NCI1), our model does not outperform the baseline, which may be a limitation of our work and requires further investigation in the future.  
experiment statistical significance: 
 Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments? Answer: [Yes] Justification: The paper reports information about the statistical significance of the experiments.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper. • The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. • It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: To see the advantage of our GRDL, we compare it with GIN. The only difference between GRDL and GIN is that GRDL uses a reference layer while GIN uses readout. We add an r ′ -layer MLP as the classifier after message-passing modules in the GIN. The following theorem gives an upper bound of GIN's generalization error:Theorem A.1 (Generalization bound of GIN). Let n = min i n i , c = ∥ Ã∥ σ , and d = max i,l d (l) i . Denote R G := c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 ( r i=1 κ (l) i ) 2 L l=1 r i=1 b (l) i κ (l) i 2/3 3 . For graphs G = {(G i , y i )} N i=1drawn i.i.d from any probability distribution over X × {1, . . . , K}, with probability at least 1 -δ, GIN network with mean readout satisfiesL γ (f ) ≤ Lγ (f ) + 3γ ln (2/δ) 2N + 8γ+24µ( r ′ i=1 κ (L+1) i ) √ R G +R ′ G ln N N where R ′ G = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 ( r i=1 κ (l) i ) 2 3C 2 2 C 1 + 3C 2 C 2 1 + C 3 1 , C 1 = r ′ i=1 b (L+1) i κ (L+1) i ,andC 2 = L l=1 r i=1 b (l) i κ (l) i 2/3 .This bound is derived using the same techniques as our GRDL bound in Theorem 3.2. To compare these two bounds, we only need to compare the following two terms:Q GRDL := 24 √ v 1 + v 2 ln N + 24γ N v 2 ln v 3 (Theorem 3.2) Q GIN := 24µ r ′ i=1 κ (L+1) i R G + R ′ G ln N (Theorem A.1)wherev 1 = 64θKµ 2 n R G , v 2 = Km d, and v 3 = 24 √ θN b D µ √ m. Our observations are as follows.• The θ in v 1 can be absorbed intoW (L) r and {D k } K k=1 . Since n ≫ K, we conclude that 64θK n is smaller than 1 in practice. Therefore, v 1 ≤ µ 2 R G . • v 2 = Km d and v 3 = 24 √ θN b D µ √ m are much smaller than R ′ G as well as R G , i.e., v 2 ≪ R ′ G and v 3 ≪ R ′ G .The reason is that R ′ G and R G involve the multiplication of terms related to c, ∥X∥ 2 2 , and κ(l) i . • In Theorem A.1, r ′ i=1 κ (L+1) iis typically larger than 3 for r ′ > 1 based on empirical observations. We also observe that r ′ i=1 κ (L+1) i may be smaller than 1 for r ′ = 1, but the linear classifier's expressive capacity is very limited and thus has large training error. Therefore, we focus on the case where r ′ > 1.Now we can conclude that Q GRDL < Q GIN . Therefore, the generalization error upper bound of GIN is larger than of GRDL, meaning our GRDL generalizes better on unseen data than GIN in the worst case. It is worth noting that these results apply to other GNNs such as GCN.We now use numerical experiments on real datasets to support our claim. The training and testing accuracy of GRDL and GIN are shown in Table 4 . We see that the training accuracy of our GRDL is close to that of GIN, but the testing accuracy of our GRDL is much higher than that of GIN. This means that GRDL and GIN have similar training errors but the former has a stronger generalization ability.Currently, we use K reference distributions for classification (one for each class). One natural approach to enhance the model's expressive power is to increase the number of reference distributions for each class. However, counterintuitively, our empirical observations suggest that having only one reference per class is optimal. In this section, we will explore and provide insights into this phenomenon. , whereD (p)k is the p-th reference in the k-th class. The prediction in Equation ( 5) is changed toEQUATION)The training algorithm is nearly the same as our GRDL with one reference per class (Algorithm 1) except for the mini-batch training loss because of the multiple references, i.e.,EQUATIONWe compare the model with P = 2 (GRDL-2) and P = 3 (GRDL-3) with our GRDL (P = 1).Table 5 shows the classification accuracy of the models on the benchmark datasets. To explain why GRDL performs better than the models with more references, we first introduce the following theorem L γ (f ) ≤ Lγ (f ) + 3γ ln (2/δ) 2N + 8γ + 24 √ v 1 + v 2 ln N + 24γ √ N v 2 ln v 3 N where v1 = 64θP 2 KRGµ 2 n , v2 = Km d, v3 = 24P √ θN bDµ √ m , RG = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 r i=1 κ (l) i 2   L l=1 r i=1 b (l) i κ (l) i 2/3   3 .This is essentially a more general version of Theorem 3.2. The following is a brief proof of this theorem.Proof. The only difference between multiple reference distributions and a single reference distribution comes from the calculation of s ij .|s ij -s ′ ij | = P p=1 MMD 2 H i , D (p) j -MMD 2 H ′ i , D (p) ′ j ≤ P p=1 MMD 2 H i , D (p) j -MMD 2 H ′ i , D (p) ′ j ≤ 4 √ θP n -1/2 ∥H i -H ′ i ∥ 2 + m -1/2 ∥D j -D ′ j ∥ 2Then with minor modifications of proof of Lemma F.14, the covering number of F is given byln N (ϵ, F, ρ) ≤ 64θP 2 KR G nϵ 2 + Kmd ln 24b D P √ θN √ mϵ .Then the theorem can be proved using the same process as the proof of Theorem 3.2.Let κ = max i,l κ(l) i , b = max i,l b (l) i κ (l) iand suppose δ is sufficiently large. The bound in Theorem B.1 can be simplified toL γ (f ) ≤ Lγ (f ) + Õ √ v 1 + γ √ N v 2 N ≤ Lγ (f ) + Õ µ b∥X∥ 2 c L (Lr) 3 2 κLr P θK/n N + γ Kmd N (26)Empirically, we observe that the training loss Lγ and the misclassification rate are nearly the same for small P and large P as shown in Figure 4 . Therefore, smaller P implies tighter generalization bound in (26). This means that one reference distribution for each class (P = 1) may be the optimal choice. Another explanation is that, the nodes embeddings of graphs in the same class can be regarded as samples drawn from a single discrete distribution, thus learning a single reference distribution is sufficient to provide high classification accuracy. On the other hand, the union of multiple distributions can be regarded as a single distribution.In practice, since the scale of θ is different from the scale of other parameters in the model, a different learning rate is used to update it. Suppose Adam [Kingma and Ba, 2014] is used to optimize the parameters, then the training of GRDL model is presented in Algorithm 1.The statistics of the datasets are reported in Table 6 . PC-3, MCF-7, and ogbg-molhiv are three large graph datasets. s i ← f w,D,θ (G i ), i ∈ B 6: ŷi ← softmax(s i ), i ∈ B 7: L = - 1 B i∈B K k=1 y ik log ŷik + λ k ′ ̸ =k ξ θ (D k , D k ′ ) 8:(g w , g D , g θ ) ← ∇ w,D,θ L 9:(w, D) ← (w, D) + α 1 • Adam ((w, D), (g w , g D ))10:θ ← θ + α 2 • Adam(θ, g θ ) 11: until convergence conditions are met 12: Output: f w,D,θIn the GNN literature, researchers typically perform 10-fold cross-validation and report the best average accuracy along with standard deviation [Xu et al., 2018 , Papp et al., 2021 , Maron et al., 2019] . But here, we adopt a different strategy used in [Vincent-Cuaz et al., 2022] . We quantify the generalization capacities of models by performing a 10-fold cross-validation with a holdout test set which is never seen during training. The validation accuracy is tracked every 5 epochs, and the model that maximizes the validation accuracy is retained for testing. This setting is more realistic than a simple 10-fold CV and can better quantify models' generalization performances [Bengio and Grandvalet, 2003] . However, the test sets for MUTAG and PTC-MR contain only 18 and 34 graphs respectively, making them too small for assessing generalization ability. Therefore, for MUTAG and PTC-MR, we use 10-fold cross-validation following [Xu et al., 2018] . Notice that kernel-based methods do not require a stopping criterion dependent on a validation set, so we report results of 10-fold nested cross-validation repeated 10 times following [Titouan et al., 2019] .Model For all the baselines, we adopt the hyper-parameters suggested in the original papers. For our methods, we use GIN [Xu et al., 2018] layers as the embedding network. Every GIN layer is an MLP of 2 layers (r = 2) with batch normalization, whose number of units is validated in {32, 64} for all datasets. The parameter λ is validated in {0.1, 1}. We validate the number of GIN layers (L) in {3, 4, 5, 6, 7, 8, 9}. For each dataset, the references' dimension (m) is validated in the minimum number of nodes (G1), average of the minimum and median number of nodes (G2), median number of nodes (G3), average of the median and maximum number of nodes (G4), and the maximum number of nodes (G5) of graphs in the dataset. The reference for each class is initialized as follows: 1) If there is a m-node graph in the dataset that belongs to the corresponding class, then the reference is initialized as node embeddings of the graph. 2) If no graph in the class has m nodes, then we perform K-Means clustering on the graphs of the class, and the m clustering center is chosen to be the initial reference.Optimization The models are trained with Adam optimizer with an initial learning rate α 1 = 10 -3 for network weights and references. The learning rate α 1 decays exponentially with a factor 0.95. Since the Gaussian kernel parameter θ is small in practice (around 10 -3 in our model), it is hard to choose a learning rate for it. Therefore, we consider π = 1 θ instead. π is initialized to 500 and its initial learning rate α 2 is validated in {0.1, 1}. The batch size for all datasets is fixed to 32.Detailed hyper-parameters setting can be found in Table 7 . To provide a comprehensive understanding, we first show the forward propagation time complexity of the three models for a single graph G = (A, X), A ∈ R n×n , X ∈ R n×d . Given that all three models employ a Graph Neural Network (GNN) for obtaining node embeddings, we denote the complexity and the number of parameters of the GNN embedding part as C 1 and N 1 , respectively. Additionally, both OT-GNN and TFGW are augmented with an MLP for classification, introducing an extra complexity C 2 and additional parameters N 2 . In the case of our GRDL, the number of references aligns with the number of classes K.Let's assume the number of references for OT-GNN is K 1 , and for TFGW, it is K 2 . Notably, TFGW and OT-GNN usually choose K 1 = 2K and K 2 = 2K. Additionally, assuming that all references D i ∈ R m×d . Since the GW distance is iteratively computed in Vincent-Cuaz et al.[2022], we denote the number of iterations for convergence as T . The time complexity and the number of parameters are outlined in Table 8 .Table 8 : Time complexity and number of parameters for GRDL, OT-GNN and TFGW.MODEL COMPLEXITY PARAMETERS GRDL C1 + O K(n 2 + mn + m 2 ) N1 + Kmd OT-GNN C1 + O K1(m + n) 3 log (m + n) + C2 N1 + K1md + N2 TFGW C1 + O T K2(m 2 n + n 2 m) + C2 N1 + K2md + N2Notice that Wasserstein distance can be approximated by sinkhorn iterations with complexity O((m + n) 2 ) per iteration [Cuturi, 2013] . But in practice, the exact calculation with O((m + n) 3 log (m + n)) complexity empirically gives better performance in terms of both precision and speed, so it is implemented in the original paper of OT-GNN [Chen et al., 2021] . Theoretically, our GRDL has lower prediction time complexity and a reduced parameter count. Table 9 shows the empirical prediction time per graph We consider two variants of GRDL. The first one is GRDL with a fixed θ = 10 -2 in the Gaussian kernel. The other is GRDL with λ = 0 in (9), which does not have the discriminativeness constraint on the references. We also include another baseline where we first use sum pooling over node embeddings and get the graph embedding vectors. The graph vectors are then used to compare with references (vectors in this case) with discrimination loss. The classification results of benchmark datasets are given in Table 11 . The original model with learnable θ and discriminativeness constraints consistently outperforms the other two.In our model, we initialized the Gaussian kernel hyper-parameter θ to 2 × 10 -3 and it was adaptively learned during training. Actually, all values between 1 × 10 -4 and 1 × 10 -1 give similar performance, as it is adaptively adjusted in the training. The initialization of the neural network parameters and the reference distributions cannot guarantee that x is close to x ′ . If θ is too large, the Gaussian kernel k(x, x ′ ) = exp (-θ∥x -x ′ ∥ 2 2 ) will be too sharp, which will lead to almost zero values. Hence, MMD will fail to effectively quantify the distance between the embeddings and reference distributions, as shown in Table 12 . θ 1 × 10 -4 1 × 10 -3 1 × 10 -2 1 × 10 -1 1 1 × 10 1 1 × 10 2 1 ×Use moderate-size message passing network We choose the training loss to be the cross-entropy loss (λ = 0) and l γ is chosen the same as the training loss. We validate the number of GIN layers L ∈ {3, 4, 5, 6} and the number of MLP layers for each GIN layer r ∈ {2, 3, 4}. L γ is set to be the validation loss and Lγ is set to be the training loss. From Figure 5 , we have the following observations• For any fixed r, if L is increasing, the empirical risk Lγ increases, and the population risk L γ either increases (r = 3 and r = 4) or first decreases then increases (r = 2). • For any fixed L, if r is increasing, the population risk L γ first decreases then increases in most of the cases, and the empirical risk Lγ decreases in most of the cases.These observations align with our bound and support our claim that moderate-size GNN should be used.Use moderate-size references We validate the reference size (m) in our experiments on real datasets, as detailed in Appendix D.3. The results in Appendix D.3 show that moderate-size references (G2, G3, G4) provide better generalization results. Here, we present the classification results for MUTAG and PROTEINS by choosing m from more fine-grained sets. For MUTAG, m is chosen from 1, 2, . . . , 30. For PROTEINS, m is chosen from 1, 3, 5, . . . , 61. The figures in Figure 6 show that our model performs the best when a moderate m is used.Regularization on references norm barely helps Consider a model regularizing the norm of references ∥D∥ 2 (GRDL-R)min L CE + λ 1 L Dis + λ 2 ∥D∥ 2 (27)The hyper-parameter λ 2 is set as 0.01. Vector and matrix norms The ℓ 2 -norm ∥ • ∥ 2 is always computed entry-wise; thus, for a matrix, it corresponds to the Frobenius norm. The metric ρ of function spaces is defined as the ℓ 2 -norm of the difference between the outputs of functions given some input X, i.e.,EQUATION)Finally, let ∥ • ∥ σ be the spectral norm and ∥ • ∥ p,q be the (p, q) matrix norm defined by ∥A∥ p,q := ∥(∥A :,1 ∥ p , . . . , ∥A :,m ∥ p )∥ q for A ∈ R d×m . Rademacher complexity Rademacher complexity is a standard complexity measure of hypothesis function space. Given dataset G and hypothesis function space F, the Rademacher complexity is defined asEQUATIONwhere σ 1 , . . . , σ N are independent Rademacher variables.Then the bound can be derived with the help of the following lemma: Lemma E.1. Given hypothesis function space F that maps a graph G ∈ X to R K and any γ > 0, define F γ := {(G, y) → l γ (f (G), y) : f ∈ F}. Then, with probability at least 1 -δ over a sampleG of size N , every f ∈ F satisfies L γ (f ) ≤ Lγ (f ) + 2R G (F γ ) + 3γ ln (2/δ) 2N .This Lemma is a standard tool in Rademacher complexity [Mohri et al., 2018] . The only problem left is to calculate the Rademacher complexity R G (F γ ).Covering number complexity bounds Direct calculation of the Rademacher complexity is often hard and the covering number is typically used to upper bound it. V is an ϵ-cover of U with respect to some metric ϱ if for all v ∈ U , there existsEQUATION)The covering number N (ϵ, U, ϱ) is defined as the least cardinality of the subset V . With covering number, the Rademacher complexity is upper bounded by the following Dudley entropy integral: Lemma E.2 (Lemma A.5 of Bartlett et al. [2017] , reformulated). Let F γ be a real-valued function class taking values in [0, γ], and assume that 0∈ F γ . Then R G (F γ ) ≤ inf α>0 4αγ √ N + 12 N γ √ N γα ln N (ϵ, F γ , ρ) dϵ .Now the only thing left is to bound N (ϵ, F γ , ρ).F Proof for Theorems, Corollaries, and LemmasWe first provide the formal definition of correct classification: This section proves some useful lemmas related to functions' lipschitz property.Lemma F.3. For any Z, Z ′ ∈ R n×d and W ∈ R d×m , ∥(Z -Z ′ )W∥ 2 ≤ ∥W ⊤ ∥ σ ∥Z -Z ′ ∥ 2Proof. First consider matrices X, Y ∈ R d×d where X, Y are positive semi-definite. Y is unitarily diagonalizable, means QΛQ -1 where Λ = diag(λ 1 , . . . , λ d ) is the diagonal matrix of eigenvalus of Y. Then we have Tr(XY) = Tr(XQΛQ -1 ) = Tr(ΛQ -1 XQ)Let P = Q -1 XQ, and λ 0 = max i λ i , we have A.3 of Bartlett et al. [2017] ). For every j and everyTr(XY) = Tr(ΛP) = k λ k p kk ≤ k λ 0 p kk = λ 0 Tr(P) = λ 0 Tr(Q -1 XQ) = λ 0 Tr(X) Take X = (Z -Z ′ ) ⊤ (Z -Z ′ ) and Y = WW ⊤ , easy to see that ∥(Z-Z ′ )W∥ 2 = Tr (Z -Z ′ ) ⊤ (Z -Z ′ )WW ⊤ ≤ λ max (WW ⊤ )Tr (Z -Z ′ ) ⊤ (Z -Z ′ ) = ∥W ⊤ ∥ σ ∥Z-Z ′ ∥ 2 Lemma F.4. If σ : R d → R d is κ-Lipschitz along ever coordinate, then it is κ-Lipschitz according to ∥ • ∥ p for any p ≥ 1. Proof. For any z, z ′ ∈ R d , ∥σ(z) -σ(z ′ )∥ p = i |σ(z) i -σ(z ′ ) i | p 1/p ≤ i κ p |z i -z ′ i | p 1/p = κ∥z -z ′ ∥ p Lemma F.5 (LemmaEQUATIONThen the proof is trivial be Lemma F.5.Lemma F.7. Cross entropy loss l(x, y) = -K k=1 y k log exp (x k ) K j=1 exp (xj ) is √ 2-Lipschitz w.r.t ∥ • ∥ 2 .Proof. Since l is differentiable, it is sufficient to find µ such that ∥∇l∥ 2 ≤ µ. Without loss of generality, supposey i = 1, then l(x, y) = -log exp (xi) K j=1 exp (xj ) . Let s = k exp x k , we have ∇l i = - l̸ =i exp (x l ) s , ∇l j = exp (x j ) s ∀j ̸ = i Therefore ∥∇l∥ 2 2 = ( l̸ =i exp (x l )) 2 + l̸ =i exp (2x l ) s 2 (a) ≤ 2 ( l̸ =i exp (x l )) 2 s 2 ≤ 2where (a) is because Cauchy-Schwarz inequality.The following lemma provides an upper bound of the covering number for the network F G .Lemma F.8 (Covering number bound of F G ). Let c = ∥ Ã∥ σ and d = max i,l d (l) i . Given an L-layer GIN message passing network F G , for any ϵ > 0 ln N (ϵ, F G , ρ) ≤ R G ϵ 2 where R G = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 κ 2 l L l=1 (τ l ) 2 3 3 and κ l = r j=1 κ (l) j , τ l = r i=1 b (l) i κ (l) i 2/3 3/2.Firstly, we introduce the core lemma used to find the covering number of compositions of multiple hypothesis function classes. Lemma F.9. Given hypothesis function classes F 1 , F 2 , . . . F k that maps input from matrix space to output in matrix space and their covering radius (ϵ 1 , ϵ 2 , . . . , ϵ k ). Assume all functions in F i is κ i -Lipschitz w.r.t. ∥ • ∥ 2 , and ln N (ϵ i , F i , ρ) ≤ g(η i ) for some function g with parameters η i (η i can be multi-valued). Then there exists ϵ-coverC of F = F k •F k-1 •• • • F 1 with ϵ = k i=1 ϵ i k j=i+1 κ j such that ln |C| ≤ k i=1 g(η i )Proof. Inductively construct covers as follows.• Let C 1 be the ϵ 1 -cover of F 1 . By our assumption,ln |C 1 | ≤ g(η 1 ) • Let C j as a ϵ j -cover of F j • • • • • F 1 . Suppose ln |C j | ≤ j i=1 g(η i ). For every f ′ j • . . . f ′ 1 ∈ C j , we construct C j+1,f ′ j ,...,f ′ 1 as an ϵ j+1 -cover of F j+1 • f ′ j • • • • • f ′ 1 . Define C j+1 := f ′ h ∈C h ,h≤j C j+1,f ′ j ,...,f ′ 1 It is clearly a cover of F j+1 • F j • • • • • F 1 By our assumption, we know that ln |C j+1,f ′ j ,...,f ′ 1 | ≤ g(η j+1 ) |C j+1,f ′ j ,...,f ′ 1 | ≤ exp (g(η j+1 )) Then, it is trivial to see |C j+1 | ≤ |C j | exp (g(η j )) ln |C j+1 | ≤ j+1 i=1 g(η i )By the inductive arguments above, we can conclude thatln |C| = ln |C k | ≤ k i=1 g(η i )Next, inductively find the cover radius ϵ for C.• It is trivial in the base case that the cover of C 1 is ϵ 1 .• Suppose for C h , the cover radius satisfiesϵ h = h i=1   ϵ i h j=i+1 κ j   For all f h+1 • f h • • • • • f 1 ∈ F h+1 • F h • • • • • F 1 , there exists f ′ h+1 • f ′ h • • • • • f ′ 1 ∈ C h+1 such that ρ(f ′ h+1 • f ′ h • • • • • f ′ 1 , f h+1 • f h • • • • • f 1 ) ≤ ρ(f ′ h+1 • f ′ h • • • • • f ′ 1 , f h+1 • f ′ h • • • • • f ′ 1 ) + ρ(f h+1 • f ′ h • • • • • f ′ 1 , f h+1 • f h • • • • • f 1 ) ≤ ϵ h+1 + κ h+1 ρ(f ′ h • • • • • f ′ 1 , f h • • • • • f 1 ) ≤ ϵ h+1 + κ h+1 ϵ h = h+1 i=1   ϵ i h+1 j=i+1 κ j   By the inductive arguments, ϵ = ϵ k = k i=1   ϵ i k j=i+1 κ j  The following matrix covering number is well-known and the detailed proof can be found in Bartlett et al. [2017] .Lemma F.10. Let conjugate exponents (p, q) and (r, s) be given with p ≤ 2, as well as positive reals (a, b, ϵ) and positive integer m. Let matrix X ∈ R n×d be given with ∥X∥ p ≤ b. Thenln N XA : A ∈ R d×m , ∥A∥ q,s ≤ a , ϵ, ∥ • ∥ 2 ≤ a 2 b 2 m 2/r ϵ 2 ln (2dm)For the composition of a hypothesis function class and a κ-Lipschitz function, we have the following lemmaLemma F.11. Suppose ψ is a κ-Lipschitz function, then ln N (ϵ, ψ • F, ρ) ≤ ln N (ϵ/κ, F, ρ) Proof. Let C denote ϵ κ -cover of F, then for any f ∈ F, there exists f ′ ∈ C such that ρ(f, f ′ ) ≤ ϵ κ . Let Z denote the input, we have ρ(ψ • f ′ , ψ • f ) = ∥ψ • f ′ (Z) -ψ • f (Z)∥ 2 ≤ κ∥f ′ (Z) -f (Zzh)∥ 2 (Lemma F.4) = κρ(f ′ , f ) ≤ ϵ Lemma F.12. Given H, H ′ ∈ R n×d and D, D ′ ∈ R m×d . Squared MMD distance with gaussian kernel k(x, y) = exp -θ∥x -y∥ 2 2 satisfies MMD 2 (H, D) -MMD 2 (H ′ , D ′ ) ≤ 4 √ θ n -1/2 ∥H -H ′ ∥ 2 + m -1/2 ∥D -D ′ ∥ 2Proof. The matrices have the formH =      h ⊤ 1 h ⊤ 2 . . . h ⊤ n      H ′ =      h ′ ⊤ 1 h ′ ⊤ 2 . . . h ′ ⊤ n      D =      d ⊤ 1 d ⊤ 2 . . . d ⊤ m      D ′ =      d ′ ⊤ 1 d ′ ⊤ 2 . . . d ′ ⊤ m     Then we haveMMD 2 (H, D) -MMD 2 (H ′ , D ′ ) ≤ 1 n 2 n i,j=1 exp -θ∥h i -h j ∥ 2 2 -exp -θ∥h ′ i -h ′ j ∥ 2 2 + 1 m 2 m i,j=1 exp -θ∥d i -d j ∥ 2 2 -exp -θ∥d ′ i -d ′ j ∥ 2 2 + 2 mn n i=1 m j=1 exp -θ∥h i -d j ∥ 2 2 -exp -θ∥h ′ i -d ′ j ∥ 2 2 (a) ≤ √ θ n 2 n i,j=1 |∥h i -h j ∥ 2 -∥h ′ i -h ′ j ∥ 2 | + √ θ m 2 m i,j=1 |∥d i -d j ∥ 2 -∥d ′ i -d ′ j ∥ 2 | + 2 √ θ mn n i=1 m j=1 |∥h i -d j ∥ 2 -∥h ′ i -d ′ j ∥ 2 | (b) ≤ √ θ n 2 n i,j=1 ∥ (h i -h ′ i ) -(h j -h ′ j ) ∥ 2 + √ θ m 2 m i,j=1 ∥ (d i -d ′ i ) -(d j -d ′ j ) ∥ 2 + 2 √ θ mn n i=1 m j=1 ∥ (h i -h ′ i ) -(d j -d ′ j ) ∥ 2 ≤ 4 √ θ n n i=1 ∥h i -h ′ i ∥ 2 + 4 √ θ m m j=1 ∥d i -d ′ i ∥ 2 (c) ≤ 4 √ θn n ∥H -H ′ ∥ 2 + 4 √ θm m ∥D -D ′ ∥ 2 =4 √ θ n -1/2 ∥H -H ′ ∥ 2 + m -1/2 ∥D -D ′ ∥ 2 .In the above derivation, (a) holds due to exp (-x 2 ) -exp (-y 2 ) ≤ |x -y| for any x, y ≥ 0, (b) holds due to the triangle inequality, and (c) holds by the Cauchy-Schwarz inequality.The covering number for a single GIN message passing layer F l with the following Lemma:Lemma F.13 (Covering number of F l ). Let c = ∥ Ã∥ σ . For any l ∈ [L] and ϵ > 0 ln N ϵ, F l , ρ ≤ c 2l τ 2 l l i=1 κ i 2 ϵ 2 ∥X∥ 2 2 ln(2 d2 ),whereκ i = j≤r κ (i) j , τ l = r i=1 b (l) i κ (l) i 2/3 3/2 .Proof. With a little abuse of notation, remove the superscript in Equation ( 19) for nowF l = {( Ã, H) → σ • • • σ ÃH W 1 • • • W r-1 W r : W i ∈ B i }whereB i := W i : ∥W ⊤ i ∥ σ ≤ κ i , ∥W i ∥ 2,1 ≤ b i . Denote F i = {Z → σ(ZW i ) : W i ∈ B i } for i ∈ [r -1], F r = {Z → ZW r : W r ∈ B r }, then F l = F r • F r-1 • • • • F 1 For any f i ∈ F i , i ∈ [r -1] with arbitrary input Z, Z ′ ∥f i (Z) -f i (Z ′ )∥ 2 = ∥σ(ZW i ) -σ(Z ′ W i )∥ 2 ≤ ∥ZW i -Z ′ W i ∥ 2 (Lemma F.4) ≤ ∥W ⊤ i ∥ σ ∥Z -Z ′ ∥ 2 (Lemma F.3) = κ i ∥Z -Z ′ ∥ 2 Similarly, for any f r ∈ F ′ r , Lemma F.3 gives ∥f r (Z) -f r (Z ′ )∥ 2 = ∥ZW r -Z ′ W r ∥ 2 ≤ κ r ∥Z -Z ′ ∥ 2Denoting Z i-1 as the input to F i (Z 0 = Z = ÃH) and using the Lipschitz conditions, we haveEQUATIONfor anyf i • f i-1 • • • • • f 1 ∈ F i • F i-1 • • • • • F 1 . So ∥Z i-1 ∥ 2 ≤ c i-1 j=1 κ j ∥H∥ 2 ≜ c i-1. By Lemma F.10 and Lemma F.11, we haveln N (ϵ i , F i , ρ) ≤ b 2 i c 2 i-1 ϵ 2 i ln(2 d2 )W is the maximum dimension of weight matrices (as previously defined in the main text). Thus by Lemma F.9, we have the covering numberln N ϵ, F l , ρ ≤ r i=1 b 2 i c 2 i-1 ϵ 2 i ln(2 d2 ) = c 2 ∥H∥ 2 2 ln(2 d2 ) r i=1 b 2 i i-1 j κ j 2 ϵ 2 i with cover radius ϵ = r i=1 ϵ i r j=i+1 κ j .Next we need to choose ϵ i to minimize the right hand side of the above inequality. Holder's inequality states that when1 p + 1 q = 1, ⟨a, b⟩ ≤ ∥a∥ p ∥b∥ q i a i b i ≤ i a p i 1/p i b q i 1/q Let α 2 i = b 2 i i-1 j κ j 2 , β i = r j=i+1 κ j . Choose p = 1 3 , q = 2 3 , i α i ϵ i 2 3 ×3 1 3 i (β i ϵ i ) 2 3 × 3 2 2 3 ≥ i (α i β i ) 2 3 i α i ϵ i 2 i β i ϵ i 2 ≥ i (α i β i ) 2 3 3    r i=1 b 2 i i-1 j=1 κ j 2 ϵ 2 i      r i=1   ϵ i r j=i+1 κ j     2 ≥ r j=1 κ 2 j r i=1 b i κ i 2 3 3 r i=1 b 2 i i-1 j=1 κ j 2 ϵ 2 i ≥ 1 ϵ 2 r j=1 κ 2 j r i=1 b i κ i 2 3 3 (33)Add the superscript back,ln N ϵ, F l , ρ ≤ ln(2 d2 ) c 2 ∥H (l-1) ∥ 2 2 ϵ 2 r j=1 κ (l) j 2   r i=1 b (l) i κ (l) i 2 3   3 H (l-1) = f (l-1) Ãf (l-2) . . . f (1) ÃX where f (k) ∈ F k for k ∈ [l -1]. By Equation ( 32), it is easy to see∥H (l-1) ∥ 2 ≤ c r j=1 κ (l-1) j ∥f (l-2) . . . f (1) ÃX ∥ 2 ≤ . . . ≤ c l-1 ∥X∥ 2 l-1 i=1   r j=1 κ (i) j   = c l-1 ∥X∥ 2 i≤l-1,j≤r κ (i) j Letting κ i = j≤r κ (i) j , τ l = r i=1 b (l) i κ (l) i 2/3 3/2, we finish the proof, i.e.,ln N ϵ, F l , ρ ≤ c 2l τ 2 l l i=1 κ i 2 ϵ 2 ∥X∥ 2 2 ln(2 d2 )With the covering number of F G , we can calculate the covering number of F in Equation ( 10). Lemma F.14 (Covering number of F). Suppose θ in the kernel (Equation (15)) is fixed. For anyϵ > 0 ln N (ϵ, F, ρ) ≤ 64θKR G nϵ 2 + Kmd ln 24b D √ θN √ mϵwhere R G is defined the same as Lemma F.8.Proof. Denote S ∈ R N ×K as the output of function F. Consider the entry (i, j) of S, by Lemma F.12, it has|s ij -s ′ ij | = MMD 2 (H i , D j ) -MMD 2 H ′ i , D ′ j ≤ 4 √ θ n -1/2 ∥H i -H ′ i ∥ 2 + m -1/2 ∥D j -D ′ j ∥ 2Then for the whole matrix S,∥S -S ′ ∥ 2 = N i=1 K j=1 |s ij -s ′ ij | 2 ≤ 4 √ θ N i=1 K j=1 n -1/2 ∥H i -H ′ i ∥ 2 + m -1/2 ∥D j -D ′ j ∥ 2 2 (a) ≤ 4 √ θ N i=1 K j=1 2 n -1 ∥H i -H ′ i ∥ 2 2 + m -1 ∥D j -D ′ j ∥ 2 2 = 4 √ 2θ (Kn -1 ∥H -H ′ ∥ 2 2 + N m -1 ∥D -D ′ ∥ 2 2 )where (a) holds due to(x + y) 2 ≤ 2(x 2 + y 2 ). If ∥H -H ′ ∥ 2 ≤ ϵ 1 and ∥D -D ′ ∥ 2 ≤ ϵ 2 , then ∥S -S ′ ∥ 2 ≤ 4 √ 2θ (Kn -1 ϵ 2 1 + N m -1 ϵ 2 2 ) ≤ ϵ by choosing ϵ 1 = √ n 8 √ Kθ ϵ and ϵ 2 = √ m 8 √ N θ ϵ. Let B D := {D ∈ R Km×d : ∥D∥ 2 ≤ b D }.It is well-known that there exists an ϵ 2 -cover obeyingN (ϵ 2 , B D , ∥ • ∥ 2 ) ≤ 3b D ϵ 2 Kmd (34)Denote the output space of F G , F as H, Z respectively, we can bound the covering number asN (ϵ, F, ρ) = N (ϵ, Z, ∥ • ∥ 2 ) ≤ N (ϵ 1 , H, ∥ • ∥ 2 ) N (ϵ 2 , B D , ∥ • ∥ 2 )It follows from Lemma F.8 and inequality (34) thatln N (ϵ, F, ρ) ≤ ln N (ϵ 1 , H, ∥ • ∥ 2 ) + ln N (ϵ 2 , B D , ∥ • ∥ 2 ) ≤ 64θKR G nϵ 2 + Kmd ln 24b D √ θN √ mϵ F.4 Proof of Lemma E.2The Dudley entropy integral bound used by Bartlett et al. [2017] is Lemma F.15 (Lemma A.5 of Bartlett et al. [2017] ). Let F be a real-valued function class taking values in [0, 1], and assume that 0 ∈ F. ThenR G (F) ≤ inf α>0 4α √ N + 12 N √ N α ln N (ϵ, F, ρ) dϵ .Lemma E.2 can be proved with simple modifications.Proof. Let F ′ = ψ • F γ where ψ(x) = 1 γ x, then F ′ is a real-valued function class taking values in [0, 1]. By Lemma F.15, it has R G (F ′ ) ≤ inf α>0 4α √ N + 12 N √ N α ln N ϵ γ , F ′ , ρ d( ϵ γ ) = inf α>0 4α √ N + 12 N γ γ √ N γα ln N ϵ γ , ψ • F γ , ρ dϵ ≤ inf α>0 4α √ N + 12 N γ γ √ N γα ln N (ϵ, F γ , ρ) dϵ (Lemma F.11)Multiplying both side by γ, it hasR G (F γ ) ≤ inf α>0 4αγ √ N + 12 N γ √ N γα ln N (ϵ, F γ , ρ) dϵ .F.5 Proof of Lemma F.8Proof. By Lemma F.13, we knowln N ϵ l , F l , ρ ≤ c 2l τ 2 l l i=1 κ i 2 ϵ 2 l ∥X∥ 2 2 ln(2 d2 )For any f l ∈ F l , given input Z, it has∥f l (Z)∥ 2 = ∥σ • • • σ ÃZ W (l) 1 • • • W (l) r-1 W (l) r ∥ 2 ≤ κ (l) r ∥σ • • • σ ÃZ W (l) 1 • • • W (l) r-1 ∥ 2 (Lemma F.3) ≤ ∥σ • • • σ ÃZ W (l) 1 • • • W (l) r-2 W (l) r-1 ∥ 2 (Lemma F.4) ≤ . . . ≤   r j=1 κ (l) j   ∥AZ∥ 2 ≤ c   r j=1 κ (l) j   ∥Z∥ 2 = cκ l ∥Z∥ 2which means all f l ∈ F l is cκ l -Lipschitz. Applying Lemma F.9, we have , β l = L j=l+1 cκ j . Choose p = 1 3 , q = 2 3 , and apply Holder's inequality in the same ways as Equation ( 33), this yieldsl α l ϵ l 2 l β l ϵ l 2 ≥ l (α l β l ) 2 3 3    L l=1 τ 2 l l j=1 cκ j 2 ϵ 2 l      L l=1   ϵ l L j=l+1 cκ j     2 ≥ c 2L L j=1 κ 2 j L l=1 (τ l ) 2 3 3 L l=1 τ 2 l l j=1 cκ j 2 ϵ 2 l ≥ 1 ϵ 2 c 2L L j=1 κ 2 j L l=1 (τ l ) 2 3 3 Thus derives the conclusion ln N (ϵ, F G , ρ) ≤ R G ϵ 2 where R G = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 κ 2 l L l=1 (τ l ) 2 3 3 .F.6 Proof of Theorem 3.2Proof. Since l γ (•, y) is µ-Lipschitz, we can bound the covering number of F γ (defined in Lemma E.1)ln N (ϵ, F γ , ρ) ≤ ln N ϵ µ , F, ρ (Lemma F.11) ≤ 64θKR G µ 2 nϵ 2 + Kmd ln 24b D µ √ θN √ mϵ (Lemma F.14) Denote v 1 = 64θKR G µ 2 n , v 2 = Km d, v 3 = 24 √ θN b D µ √ m, then by Lemma E.2, we can bound the Rademacher complexityR G (F γ ) ≤ inf α>0 4αγ √ N + 12 N γ √ N γα v 1 ϵ 2 + v 2 ln v 3 ϵ dϵ (a) ≤ inf α>0 4αγ √ N + 12 N γ √ N γα v 1 + v 2 ϵ 2 + v 2 ln v 3 dϵ (b) ≤ inf α>0 4αγ √ N + 12 N γ √ N γα √ v 1 + v 2 ϵ + v 2 ln v 3 dϵ = inf α>0 4αγ √ N + 12 N √ v 1 + v 2 ln √ N α + γ v 2 ln v 3 √ N -α (c) ≤ 4γ N + 12 √ v 1 + v 2 ln N N + 12γ(N -1) √ v 2 ln v 3 N √ N ≤ 4γ N + 12 √ v 1 + v 2 ln N N + 12γ √ v 2 ln v 3 √ Nwhere (a) holds due to ln 1 x ≤ 1 x 2 , (b) holds due to √ x + y ≤ √ x + √ y.L γ (f ) ≤ Lγ (f ) + 8γ + 24 √ v 1 + v 2 ln N + 24γ √ N v 2 ln v 3 N + 3γ ln (2/δ) 2N This finished the proof.EQUATIONLζ (f ) ≤ N -1 i 1[f (G i ) yi ≤ ζ + arg max j̸ =yi f (G i ) j ].Proof. Choose the loss l γ (•, y) asl γ (•, y) = r ζ (-M(•, y))where M(v, y) := v y -max i̸ =y v i is the margin operator andr ζ (t) :=    0 t < -ζ, 1 + t/ζ t ∈ [-ζ, 0], 1 t > 0.is called the ramp loss. The population ramp risk is defined asL ζ (f ) := E G∼X [r ζ (-M(f (G), y)))].Given the graph dataset G sampled from X , the empirical ramp risk is Lζ (f ) :=N -1 N i=1 r ζ (-M(f (G i ), y i )). It is clear that 1[arg max j f (G) j ̸ = y] ≤ r ζ (-M(f (G), y)), so Pr G∼X [arg max j f (G) j ̸ = y] = E G∼X 1(arg max j f (G) j ̸ = y) ≤ L ζ (f )It is easy to see that γ = 1 in this case. Also by Lemma F.6, µ = 2 ζ . Then it is trivial to get the bound by Theorem 3.2 with simple substitutionPr G∼X [arg max j f (G) j ̸ = y] ≤ L ζ (f ) ≤ Lζ (f ) + 3 ln (2/δ) 2N + 8 + 24 √ v 1 + v 2 ln N + 24 √ N v 2 ln v 3 N where v 1 = 256θKR G nζ 2 , v 2 = Km d, v 3 = 48b D √ θN √ mζ , R G = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 r i=1 κ (l) i 2 L l=1 r i=1 b (l) i κ (l) i 2/3 3 . Also, r ζ (-M(f (G i ), y i )) ≤ 1[f (G i ) yi ≤ ζ + arg max j̸ =yi f (G i ) j ], so we have Lζ (f ) = N -1 N i=1 r ζ (-M(f (G i ), y i )) ≤ N -1 N i=1 1[f (G i ) yi ≤ ζ + arg max j̸ =yi f (G i ) j ]Corollary F.17 (Generalization bound with cross-entropy loss). Suppose l γ (•, y) is the cross-entropy loss L CE (7). Let n be the minimum number of nodes for graphs {G i } N i=1 , θ be the hyper-parameter in gaussian kernel (Equation (15)), c = ∥ Ã∥ σ . For graphs G = {(G i , y i )} N i=1 drawn i.i.d from any probability distribution over X × {1, . . . , K} and references {D k } K k=1 , D k ∈ R m×d , with probability at least 1 -δ, every network f ∈ F under Assumption 3.1 satisfyL γ (f ) ≤ Lγ (f ) + 3γ ln (2/δ) 2N + 8γ + 24 √ v 1 + v 2 ln N + 24γ √ N v 2 ln v 3 N where v 1 = 128θKR G n , v 2 = Km d, v 3 = 24 √ 2θN b D √ m , R G = c 2L ∥X∥ 2 2 ln(2 d2 ) L l=1 r i=1 κ (l) i 2 L l=1 r i=1 b (l) i κ (l) i 2/3 3 . Proof. According to Lemma F.7, µ = √ 2.Then the proof is trivial by substituting it into Theorem 3.2.F.8 Adjacency matrix spectral norm Lemma F.18. Let G = (V, E) be an undirected graph with adjacency matrix A ∈ R n×n , d G be the maximum degree of G. Then, the adjacency matrix satisfies∥A∥ σ ≤ d G Proof.Based on the definition of the spectral norm, we have∥A∥ σ (a) = max ∥x∥2=1 x ⊤ Ax = max ∥x∥2=1 (i,j)∈E x i x j ≤ max ∥x∥2=1 (i,j)∈E 1 2 (x 2 i + x 2 j ) = d G i∈V x 2 i = d Gwhere (a) is because A is a real symmetric matrix.Lemma F.19. For any matrix X ∈ R m×n , ∥X∥ σ ≤ ∥X∥ 2 ≤ √ r∥X∥ σ where r = rank(X)Proof. By the definition of spectral norm,∥X∥ σ = λ max (X ⊤ X)where λ max denotes the largest eigenvalue. Since X ⊤ X is a positive semi-definite real symmetric matrix, it must has n real eigenvalues that can be ordered asλ 1 ≥ λ 2 ≥ • • • ≥ λ r > λ r+1 • • • = λ n = 0.Then it has∥X∥ σ = λ 1 ≤ n i=1 λ i = tr(X ⊤ X) = ∥X∥ 2 ≤ rλ 1 = √ r∥X∥ σ . Lemma F.20. Let G = (V, E) be a graph with adjacency matrix A ∈ R n×n . Assume |E| > 0, then c = ∥ Ã∥ σ = ∥A + I∥ σ > 1.Proof. By Lemma F.19,∥ Ã∥ σ ≥ 1 √ r ∥ Ã∥ 2 = 1 √ r ∥A + I∥ 2 ≥ 1 √ r ∥A∥ 2 + 1 √ r ∥I∥ 2 (a) > 1 √ r ∥I∥ 2 = n r ≥ 1where (a) is because |E| > 1.F.9 Generalization of MMD (to be used in Section F.10)Before stating the lemma, we first give alternative definitions of MMD. Let P be a continuous probability distribution of some random variable Z taking values from space Z. Then, the kernel mean embedding of P associated with the continuous, bounded, and positive-definite kernel function kEQUATIONwhich is an element in the Reproducing Kernel Hilbert Space (RKHS) H associated with kernel k. In many practical situations, it is unrealistic to assume access to the true distribution P. Instead, we only have access to samples P = {z i } n i=1 from P. We can approximate (35) by the empirical kernel mean embeddingEQUATIONFor another continuous distribution Q with samples Q = {z ′ i } m i=1 , the MMD between the two probability distribution is defined asMMD(P, Q) = ∥µ P -µ Q ∥ H ,and the empirical MMD is MMD(P, Q) = ∥μ P -μQ ∥ H .Denote d := MMD(P, Q) and d := MMD(P, Q), we have the follow Lemma:Lemma F.21. With probability at least 1 -δ we have|d -d| ≤ 1 √ m + 1 √ n 2 + 2 log 2 δProof. We use the following Lemma:Lemma F.22 (Theorem 7 of Gretton et al. [2012a] , reformulated). Assume 0 ≤ k(x, y) ≤ K. Then with probability at least 1 -2 exp -ε 2 mn 2K(m+n)d -d ≤ 2 1 √ m + 1 √ n + ε Let δ = 2 exp -ε 2 mn 2K(m+n) , then it has ε = 1 m + 1 n 2 log 2 δ .Therefore, with probability at least 1-δ d -d ≤ 2 1 √ m + 1 √ n + 1 m + 1 n 2 log 2 δ ≤ 1 √ m + 1 √ n 2 + 2 log 2 δ F.10 Proof of Lemma F.2Proof. For simplicity, let k ′ := arg min j̸ =k MMD(H i , D j ). Since H i and D k are finite samples from P k , and D k ′ is sampled from P k ′ , by Lemma F.21, we have|MMD(P k , P k ) -MMD(H i , D k )| ≤ ∆ 1 = 1 √ m + 1 √ n 2 + 2 log 2 δ ′ w.p. (1 -δ ′ ), |MMD(P k , P k ′ ) -MMD(H i , D k ′ )| ≤ ∆ 2 = 1 √ m + 1 √ n 2 + 2 log 2 δ ′ w.p. (1 -δ ′ ).Therefore, with probability at least 1 -2δ ′ (union bound), we haveMMD(H i , D k )-MMD(P k , P k ) = MMD(H i , D k )-0 ≤ ∆ 1 and MMD(P k , P k ′ )-MMD(H i , D k ′ ) ≤ ∆ 2 .It follows thatMMD(H i , D k ) -MMD(H i , D k ′ ) ≤ -MMD(P k , P k ′ ) + ∆ 1 + ∆ 2 .By Definition F.1, to ensure correct classification, we can let-MMD(P k , P k ′ ) + ∆ 1 + ∆ 2 < 0.This meansEQUATIONTherefore, if (37) holds, the classification is correct with probability at least 1 -2δ ′ . Letting δ = 2δ ′ , we finish the proof.Let H ∈ R i ni×d be the output of the message passing layers. Then the mean readout is a matrix multiplicationZ = MH, M =      1 n1. . . N (ϵ, F ′ , ρ) ≤ c 2L ∥X∥ 2 2 ln(2 d2 ) ϵ 2 A where A = L l=1 ( r i=1 κ (l) i ) 2 ( r ′ i=1 κ (L+1) i ) 2 L l=1 r i=1 b (l) i κ (l) i 2/3 + r ′ i=1 b (L+1) i κ (L+1) i 2/3 3 .Then by Lemma E.2, the generalization bound can be easily derived by taking α = 1 √ N .• The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide a short proof sketch to provide intuition. • Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material. • Theorems and Lemmas that the proof relies upon should be properly referenced. 


Decision: The paper proposes an architecture for graph-level prediction based on node-level representation computing by a GNN. Instead of using a pooling method, the paper derives a reference distribution learning method for graph-learning tasks. That is, the authors' architecture measures the dissimilarity between distributions of the input graph and the references. The authors propose a generalization analysis based on standard tools and conducted an empirical study showing somewhat promising results, mostly on smaller datasets from the TUDataset collection.

The reviewers liked the novel idea of using distribution comparison for graph-level prediction. Some highlighted the clear motivation and good exposition and the performed extensive ablation studies on the proposed method. Reviewer anD1 argued that the method shares similarities with mean pooling and max pooling; the time complexity is too high, and some GNN pooling baselines are missing in the experiments. The authors partially addressed some of these concerns in the rebuttal.

This is a borderline paper. A majority of reviewers favor its acceptance due to its novelty. I follow this assessment.

Reviews: review: paper_summary: paper_summary: This paper proposes to make graph-level predictions via distribution comparison between node-level representations and discrete reference distributions. The authors claim that their proposed method avoids the requirements of graph pooling for graph-level tasks. and reduce the risk of information loss. Theoretical and empirical justification results are provided.
strength_weakness: strength_weakness: - Equation in L138-139 (a missing equation index?) also requires summation over nodes. The description of "avoid graph pooling operation" seems to be overclaimed. I suggest that the authors reconsider it.
- An ablation study on the discrimination loss and the usage of node-level representation is required. I think the main improvement in performance may be attributed to the discrimination loss. It helps the model learn distant representations for graphs of different labels in the feature space.  The performance of the proposed method should be compared with and without the discrimination loss. Besides, the authors should also implement a baseline with discrimination loss where node representations are first sum-pooling and then compared with the reference distribution.
- Minor: Latest baseline models are required for a comprehensive empirical comparison, such as graph pooling including SEP[1], GMT[2], and CoCN[3] and graph transformers including Exphormer[4], GRIT[5], and MPNN-VN[6]. Considering the time limitation and the similarity of the datasets, a comparison on part of the datasets will be sufficient.

[1] Structural entropy guided graph hierarchical pooling. ICML'22

[2] Accurate Learning of Graph Representations with Graph Multiset Pooling. ICLR'22

[3] All in a row: Compressed convolution networks for graphs. ICML'23

[4] Exphormer: Sparse Transformers for Graphs. ICML'23

[5] Graph Inductive Biases in Transformers without Message Passing. ICML'23

[6] On the Connection Between MPNN and Graph Transformer. ICML'23
questions: questions: Please refer to the weaknesses.
limitations: limitations: Yes. The authors describe the limited performances of the proposed method on certain datasets.
paper_summary: paper_summary: The paper introduces a novel algorithm called GRDL for graph classification. GRDL  treats each graph’s latent node embeddings given by GNN layers as a discrete distribution and directly classify distributions without global pooling. The authors derived generalization error bounds for GRDL and verified them numerically. The experiments on many graph datasets show the superiority of GRDL. GRDL is 10 times faster than leading competitors.
strength_weakness: strength_weakness: There is no major weakness found. Please refer to my questions in the next section.
questions: questions: * Besides GIN, there are other GNN models such as GCN and GAT. Why did the authors consider GIN only in the theoretical analysis and experiments?
* I think the parameter $\theta$ of the Gaussian kernel can be absorbed into the neural network parameter. I suggest the authors discuss the necessity of setting or optimizing $\theta$ separately. 
* In Theorem 3.2, the bound is related to $K$, the number of classes. But, in many literature of generalization analysis, the bound is not related to the number of classes. Could the author explain the difference?
limitations: limitations: The authors discussed the limitations in Section 6.
paper_summary: paper_summary: Instead of compressing node embedding matrix into a graph-level vector, this paper proposes a reference distribution learning method GRDL especially designed for graph classification task. GRDL achieves graph classification by measuring the dissimilarity between distributions of the input graph and the references. Theoretical analysis of the generalization error bound offers guidance of hyperparameter tuning.
strength_weakness: strength_weakness: - The idea is borrowed from domain transfer learning but lacks adaption to graph classification task. In graph classification task, graphs in each dataset “are drawn i.i.d” as the authors mentioned at line 205. This means that although graphs have different classes, they belong to an identical distribution, which is contradictory to the hypothesis of this paper that different classes belong to different reference distributions.
- In my opinion, this paper didn’t avoid information compression compared to mean pooling and max pooling. MMD unifies the source and target inputs as a vertex in the same RKHS and then calculate the distance between them. I think the mapping from embedding matrix into a vertex can be regarded as information compression, and the calculation of mean discrepancy is similar to mean pooling.
- The time complexity is too high, the datasets used for large-scale experiments are unreasonable, and the competitors used in efficiency experiment is unreasonable. The time complexity is $O(N^2)$, which is the same as the highly criticised time-consuming node clustering pooling methods (such as DiffPool[1]), while node dropping methods (such as SAGP[2]) only requires a time complexity of $O(E)$. The complexity restricts its scalability to larger-scale graphs. Three so-called large-scale datasets are not really large-scale. They are only large in the number of graphs but not in the number of nodes in each single graph. The authors should add experiments on synthetic large-scale datasets. Besides, in time cost comparison, the authors didn’t choose competitors with SOTA efficiency but used two time consuming methods, which is quite unreasonable.
- The theoretical contribution of generalization error bound is limited. Firstly, the result offers fuzzy insights into the choose of hyperparameters by giving “moderate-size message passing GIN”, and “moderate-size references”. How to define “moderate”? At line 240, the authors analysis that “a network with a smaller $L4 and $r4 may guarantee a tighter bound on the population risk compared to a larger one. Therefore, a promising strategy is to use a moderate-size message passing GNN”. What’s the causality between “smaller” and “moderate-size”? That’s quite confusing. Experimental results in Appendix show that hyperparameters are still choosed on trial and error. Besides, the generalization ability comparison with GIN with mean readout is meaningless because the competitor is not the SOTA.
- Important baselines are missing in the experiments. Please add comparison results with recent graph pooling methods such as ASAP[3], MinCutPool[4], StructPool[5], MuchPool[6], TAP[7], Wit-TopoPool[8], and MSGNN[9].

References

[1]Ying Z, You J, Morris C, et al. Hierarchical graph representation learning with differentiable pooling[J]. Advances in neural information processing systems, 2018, 31.

[2]Lee J, Lee I, Kang J. Self-attention graph pooling[C]//International conference on machine learning. PMLR, 2019: 3734-3743.

[3]Ranjan E, Sanyal S, Talukdar P. Asap: Adaptive structure aware pooling for learning hierarchical graph representations[C]//Proceedings of the AAAI conference on artificial intelligence. 2020, 34(04): 5470-5477.

[4]Bianchi F M, Grattarola D, Alippi C. Spectral clustering with graph neural networks for graph pooling[C]//International conference on machine learning. PMLR, 2020: 874-883.

[5]Yuan H, Ji S. Structpool: Structured graph pooling via conditional random fields[C]//Proceedings of the 8th International Conference on Learning Representations. 2020.

[6]Du J, Wang S, Miao H, et al. Multi-Channel Pooling Graph Neural Networks[C]//IJCAI. 2021: 1442-1448.

[7]Gao H, Liu Y, Ji S. Topology-aware graph pooling networks[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021, 43(12): 4512-4518.

[8]Chen Y, Gel Y R. Topological pooling on graphs[C]//Proceedings of the AAAI Conference on Artificial Intelligence. 2023, 37(6): 7096-7103.

[9]Lv Y, Tian Z, Xie Z, et al. Multi-scale Graph Pooling Approach with Adaptive Key Subgraph for Graph Representations[C]//Proceedings of the 32nd ACM International Conference on Information and Knowledge Management. 2023: 1736-1745.
questions: questions: Please refer to the weaknesses.
limitations: limitations: N/A

Comments: comment: Dear Reviewer anD1,

We appreciate your comments and suggestions. Did our response address your concerns? We are keen to receive your feedback and provide further explanation if necessary.

Sincerely,

Authors
Thank you so much for your feedback and support. The suggestions from you and the other two reviewers have helped us improve the quality of our paper. We will update the manuscript according to the rebuttal.
I appreciate the detailed response which has addressed my concerns. The new ablation results (W2) further validate the proposed method and the selected comparison results seem promising. Please make sure to clarify the difference between the pre-pooling methods and your discrepancy measuring and summation strategy, and update your manuscript based on the rebuttal.
We sincerely appreciate your recognition of our work.
The rebuttal and the additional experimental results have adequately addressed my concerns. Therefore, I raise my confidence level to 5.
Sorry for the confusion. All additional experiments can be found in the global rebuttal. We did not update the PDF. Thanks for pointing it out.
The authors referred to a PDF file in their global rebuttal, but I am unable to locate it. Could the authors please confirm whether this file has been correctly uploaded?