A PIECE-WISE POLYNOMIAL FILTERING APPROACH
Year: 2021

Abstract:
Graph Neural Networks (GNNs) exploit signals from node features and the input graph topology to improve node classification task performance. However, these models tend to perform poorly on heterophilic graphs, where connected nodes have different labels. Recently proposed GNNs work across graphs having varying levels of homophily. Among these, models relying on polynomial graph filters have shown promise. We observe that solutions to these polynomial graph filter models are also solutions to an overdetermined system of equations. It suggests that in some instances, the model needs to learn a reasonably high order polynomial. On investigation, we find the proposed models ineffective at learning such polynomials due to their designs. To mitigate this issue, we perform an eigendecomposition of the graph and propose to learn multiple adaptive polynomial filters acting on different subsets of the spectrum. We theoretically and empirically show that our proposed model learns a better filter, thereby improving classification accuracy. We study various aspects of our proposed model including, dependency on the number of eigencomponents utilized, latent polynomial filters learned, and performance of the individual polynomials on the node classification task. We further show that our model is scalable by evaluating over large graphs. Our model achieves performance gains of up to 10% over the state-of-the-art models and outperforms existing polynomial filter-based approaches in general. Our anonymized code is available at https://tinyurl.com/PPGNN.

heading: 1 INTRODUCTION
text: In this work, we are interested in the problem of classifying nodes in a graph. In this problem, a graph and features for all the nodes in the graph are made available. We are also given labels for a subset of the nodes in the graph. The learning model utilizes this information to predict labels for the remaining nodes. Graph Neural Networks (GNNs) perform well on such problems (Kipf & Welling, 2017). GNNs predict a node’s label by aggregating information from its neighbours. Thus, the performance of GNN is dependent on how well the given graph correlates with the node labels. Characterizing the correlation between the graph and node labels is an active area of research. Several metrics have been proposed including edge homophily (Abu-El-Haija et al., 2019a; Zhu et al., 2020), node homophily (Pei et al., 2020), class homophily (Lim et al., 2021). All these metrics show that GNNs perform well when the graphs and node labels are positively correlated. For example, in the simplest case, GNNs work well when the node and its neighbours share similar labels. However, the performance can be poor if this criterion is not satisfied.
Several proposed GNN models attempt to be robust to the underlying correlation between graphs and labels. Some approaches modify the aggregation mechanism (Pei et al., 2020; Zhu et al., 2020; Kim & Oh, 2021), while other approaches propose to estimate and leverage the label-label compatibility matrix as a prior (Zhu et al., 2021). More recent approaches have tackled this problem from a filter learning perspective (Bo et al., 2021; Chien et al., 2021). These methods learn a filter function that operates on the eigenvalues of the graph. With eigenvalues having frequency interpretations (Shuman et al., 2013), the filter function selectively accentuates and suppresses various frequencies as required by the task. This process enables the model to learn better embeddings, which translates to improved performance. Our interest lies in the area of modelling GNNs that can learn an effective filter.
Chien et al. (2021) proposed to learn a polynomial filter. This model achieves good performance gains on several datasets with varying correlations between graphs and labels. However, we observe that learning a polynomial of a degree lower than the number of nodes can be seen as solving an over-determined system of equations. Our observation suggests needing a higher-order polynomial to model richer or more complex frequency responses for certain datasets. The design of the model proposed in Chien et al. (2021) makes it ineffective in learning a higher-order polynomial. It happens because the model attempts to overcome over-smoothing by giving smaller coefficient values to larger powers. In this paper, we address this issue and make the following contributions:
• Inspired by Chien et al. (2021), we propose to learn a polynomial filter; however, we do so by learning a sum of polynomials over different subsets of eigenvalues. Such modelling enables learning higher-order polynomials with several low order polynomials acting over the different subsets.
• Such modelling, however, requires an eigendecomposition of the graph, which can be expensive. We leverage the model in Chien et al. (2021) and the fact that there exist efficient algorithms for computing top and bottom eigen components to propose a practically efficient model.
• We present theoretical analysis that suggests that our sum of polynomials approach can approximate a latent optimal filter better than a single polynomial. We also show that the space of learnable filters and graphs using our approach is larger than what is possible with Chien et al. (2021).
• Our experimental results show that our model performs better than state-of-the-art methods on several datasets, giving up to 10% gains. These gains indicate that modelling the filter as a sum of polynomials is needed to approximate the latent filter better. We also answer other research questions through ablative studies, including the importance of the top and bottom eigenvalues, the number of top and bottom eigenvalues needed, the number of polynomials needed, and the trade-offs involved.
The rest of the paper is organized as follows: We present related work in Section 2, problem setup and motivation for the approach in Section 3, details of the proposed approach with complexity and theoretical analysis in Section 4 and finally experimental results and empirical analysis in Section 5.

heading: 2 RELATED WORKS
text: In the recent times, Graph Neural Networks (GNNs) have become an increasingly popular method for semi-supervised classification with graphs. Bruna et al. (2014) set the stage for early GNN models, which was then followed by various modifications (Defferrard et al., 2016; Kipf & Welling, 2017; Hamilton et al., 2017; Veličković et al., 2018). Incorporating random walk information (AbuEl-Haija et al., 2019b;a; Li et al., 2018) gave further improvements in these models, but they still suffer from over smoothing. Several proposed models attempt to overcome this problem (Klicpera et al., 2019; Lukovnikov & Fischer, 2021; Chamberlain et al., 2021; Yang et al., 2021).
Another line of research explored the question of which graphs worked well with GNNs. The critical understanding was that the performance of GNN was dependent on the correlation of the graphs with the node labels. Several approaches (Abu-El-Haija et al., 2019a; Zhu et al., 2020; Wang & Derr, 2021) considered edge homophily and proposed a robust GNN model by aggregating information from several higher-order hops. Kim & Oh (2021) also considered edge homophily and mitigated the issue by learning robust attention models. Pei et al. (2020) talks about node homophily and proposes to aggregate information from neighbours in the graph and neighbours inferred from the latent space. Zhu et al. (2021) proposes to estimate label-label compatibility matrix and uses it as a prior to update posterior belief on the labels. However, these approaches still rely on the given graphs for information like the hop neighbours etc., which might be poor approximations of the desired neighbours with similar labels.
Some of the recent approaches focused on learning filter functions that operate on the eigenvalues of the graph. Learning these filter functions can be viewed as directly adapting the graph for the desired task. Bo et al. (2021) models the filter function as an attention mechanism on the edges, which learns the difference in the proportion of low-pass and high-pass frequency signals. Chien et al.
(2021) proposes a polynomial filter on the eigenvalues that directly adapts the graph for the desired task. Zheng et al. (2021) decompose the graph into low-pass and high-pass frequencies, and define a framelet based convolutional model. Our work is closely related to these lines of exploration. In this work, we propose to learn a filter function as a sum of polynomials over different subsets of the eigenvalues, enabling design of effective filters to model task-specific complex frequency responses with compute trade-offs.

heading: 3 PROBLEM SETUP AND MOTIVATION
text: We focus on the problem of semi-supervised node classification on a simple graph G = (V, E), where V is the set of vertices and E is the set of edges. Let A ∈ {0, 1}n×n be the adjacency matrix associated with G, where n = |V| is the number of nodes. Let Y be the set of all possible class labels. Let X ∈ Rn×d be the d-dimensional feature matrix for all the nodes in the graph. Given a training set of nodes D ⊂ V whose labels are known, along with A and X, our goal is to predict the labels of the remaining nodes. Let AI = A + I where I is the identity matrix. Let DAI be the degree matrix of AI and Ã = D −1/2 AI AID −1/2 AI
. Let Ã = UΛUT be the eigendecomposition of Ã. Then, the spectral convolution of X on the graph A can be defined via the reference operator Ã and a filter function h operating on the eigenvalues, in the Fourier domain (Tremblay et al., 2017; Chien et al., 2021) as,
Z = UH(Λ)UTX = Udiag(Vα)UTX = k∑ j=1 αjÃ jX (1)
where H(Λ) = diag(h(Λ)). In this process, the filter function is essentially adapting the graph for the desired task at hand. The second equality follows from using a polynomial filter, H(λ) = ∑k i=1 αiλ
i where αi’s are coefficients of the polynomial, k is the order of the polynomial and λ is any eigenvalue from Λ. Note that H(Λ) = diag(Vα) where V ∈ Rn×k is a Vandermonde matrix constructed using eigenvalues from Λ and α ∈ Rk×1 is the vector consisting of αj’s. The last equality follows from the eigendecomposition properties, enabling avoidance of eigendecomposition. It is well-known that polynomial filters can approximate any graph filter (Shuman et al., 2013; Tremblay et al., 2017). However, several practical challenges arise in learning a good polynomial filter. We briefly discuss them below.
Let h∗ be the optimal filter for some given task. Learning a polynomial filter involves the approximation: h∗ ≈ Vα. Notice that this system is over-determined since the number of unknowns (k, the size of α) is less than the number of equations (n, the number of nodes). To obtain a consistent solution, we may have to work with higher-order polynomials. Working with higher powers poses several practical challenges, and the following observations are in order. First, node features computed via ÃjX becomes indistinguishable (see Figure 1a and A.2.1). Since the node features are indistinguishable, the importance of the coefficient associated with ÃjX for the task at hand reduces significantly. These observations are in line with Theorem 4.2 in Chien et al. (2021), where it states that as over-smoothing happens, the corresponding coefficients of the polynomial go to zero. We observe that with insignificant contributions from higher-order terms, the performance (aka test accuracy) does not improve with the increase in the order of the polynomial (See Figure 1b, 1c and A.2.2). Next, one can think of addressing this issue with ÃjX by directly working with full eigendecomposition. However, this approach is prohibitively expensive for large graphs. Even if we were to do this, λi’s are in the range of [−1, 1] and thus λji will diminish with higher powers. Therefore, it will continue to be less effective in learning coefficients of a higher power. Motivated by the above observations, we propose a novel piece-wise polynomial filter learning approach.

heading: 4 PROPOSED APPROACH
text: The central goal of our approach is to get improved task-specific performance by learning effective task-specific graph filters. We propose to learn the polynomial filter as a sum of polynomials operating on different eigenvalue intervals, taking task-specific requirements and practical considerations into account. We show that our proposed filter design can approximate the latent optimal graph
filter better than a single polynomial, and the resultant class of learnable filters/graphs is richer. We briefly discuss the training and computational complexity of the proposed model.

heading: 4.1 PIECE-WISE POLYNOMIAL/SPLINE GRAPH FILTER FUNCTION
text: We model the filtering function h(λ) as a piece-wise polynomial or spline function where each polynomial is of a lower degree (e.g., a cubic polynomial). We partition the eigen spectrum in [−1, 1] (or [0, 2] as needed) into contiguous intervals and approximate the desired frequency response by fitting a low degree polynomial in each interval. This process helps us to learn a more complex shaped frequency response as needed for the task. Let S = {σ1, σ2, . . . , σm} denote a partition set with m contiguous intervals and hi,ki(λ; γi) denote a ki-degree polynomial filter function of the interval σi with polynomial coefficients γi. We define PP-GNN filter function as:
h(λ) = ∑ σi∈S hi,ki(λ; γi) (2)
and learn a smooth filter function by imposing additional constraints to maintain continuity between polynomials of contiguous intervals at different endpoints (aka knots). This class of filter functions is rich, and its complexity is controlled by choosing intervals (i.e., endpoints and number of partitions) and polynomial degrees. Given the filter function, we compute the node embedding matrix as:
Z = m∑ i=1 UiHi(γi)U T i Z0(X; Θ) (3)
where Ui is a matrix with eigenvectors corresponding to eigenvalues that lie in σi, Hi(γi) is the diagonal matrix with diagonals containing the hi evaluated at the eigenvalues and Z0(X; Θ) is an MLP network with parameters Θ.

heading: 4.2 PRACTICAL AND IMPLEMENTATION CONSIDERATIONS
text: The filter function (3) requires computing eigendecomposition of Ã and is expensive, therefore, not scalable for very large graphs. Thus, having finer control to learn complex frequency responses comes with a high computational cost. We address this problem by devising a graph filter function that enables a trade-off between computational cost and frequency response control. First, we make the following observations:
1. Top-k and bottom-k eigenvalues represent k low and high frequencies. There are efficient algorithms and off-the-shelf library packages available to get top-k and bottom-k eigenvalues and eigenvectors of sparse matrices.
2. GPR-GNN method also learns a graph filter but operates on the entire eigen spectrum by sharing the filter coefficients across the spectrum. Therefore, it is a special case of our proposed model (3). Since GPR-GNN learns a global polynomial by having shared coefficients across the spectrum (like learning a global polynomial to fit a function), it is not flexible enough to learn complex frequency responses, as needed in several practical applications. Note that increasing the degree of the global polynomial does not help beyond some high degree specific to tasks. We demonstrate this through extensive experiments on benchmark datasets in the experiment and appendix sections. One important advantage of GPR-GNN is that it does not require computation of eigendecomposition; therefore, efficient.
3. Many recent works, including GPR-GNN, investigated the problem of designing robust graph neural networks that work well across homophilic and heterophilic graphs. They found that graph filters that amplify or attenuate low and high-frequency components of signals (i.e., low-pass and high-pass filters) are critical to improving performance on several benchmark datasets.
Efficient Variant. Using the above observations, we propose an efficient variant of (2) as: h̃(λ) = ηl ∑ σi∈Sl h (l) i (λ; γ (l) i ) + ηh ∑ σi∈Sh h (h) i (λ; γ (h) i ) + ηgprhgpr(λ; γ) (4)
where Sl consists of partitions over low frequency components, Sh consists of partitions over high frequency components, the first and second terms fit piece-wise polynomials1 in low and high frequency regions, as indicated through superscripts. However, we do not want to lose any useful information from other frequencies in the central region, yet maintain efficiency. We achieve this by adding the GPR-GNN filter function, hgpr(λ; γ), which is computationally efficient. Since hgpr(λ; γ) is a special case of (2) and the terms in (4) are additive, it is easy to see that (4) is same as (2) with a modified set of polynomial coefficients. In practice, the required low and high-frequency components are computed based on affordable computational cost. Furthermore, we can also control the contributions from each term by setting or optimizing over hyperparameters, ηl, ηh and ηgpr. Thus, the proposed model offers richer capability and flexibility to learn complex frequency response and balance computation costs over GPR-GNN.
Model Training. Like GPR-GNN, we apply SOFTMAX activation function on (3) and use the standard cross-entropy loss function to learn the sets of polynomial coefficients (γ) and classifier model parameters (Θ) using labeled data. To ensure smoothness of the learned filter functions, we add a regularization term that penalizes squared differences between the function values of polynomials of contiguous intervals at each other’s interval end-points. More details can be found in the appendix (A.4.1).
Computational Cost. There is some pre-training cost of computing the eigendecomposition for top and bottom k eigenvalues. Most algorithms for this task utilize Lanczos’ iteration, convergence bounds of which depends on the input matrix’ spectrum (Saad, 1980; LI, 2010), which although have superlinear convergence, but are observed to be efficient in practice. We compute node embeddings afresh whenever the model parameters are updated. This computation involves matrix multiplication with eigenmatrices incurring an additional cost (over GPR-GNN) of O(nkL) where k and L denote the number of selected low/high eigenvalues and classes, respectively.

heading: 4.3 ANALYSIS
text: We first give a simple theoretical justification for using multiple polynomials. We upper-bound the approximation error achieved by using multiple polynomials by the error that the single polynomial parameterization achieves. We relegate the proofs of all the theorems of this section to the appendix.
Theorem 4.1. For any frequency response h∗, and an integer K ∈ N, let h̃ := h + hf , with hf having a continuous support over a subset of the spectrum, σf . Assume that h and hf are parameterized by independent K and K ′-order polynomials, p and pf , respectively, with K ′ ≤ K. Then there exists h̃, such that min ‖h̃−h∗‖2 ≤ min ‖h−h∗‖2, where the minimum is taken over the polynomial parameterizations. Moreover, for multiple polynomial adaptive filters hf1 , hf2 , ..., hfm
1For brevity, we dropped the polynomial degree dependency.
parameterized by independentK ′-degree polynomials withK ′ ≤ K but having disjoint, contiguous supports, the same inequality holds for h̃ = h+ ∑m i=1 hfi .
For a detailed proof please refer to A.4.3. We constructed waveforms of arbitrary complexity, approximated by first fitting a single globally-defined polynomial and then adding locally defined polynomials on top of it. We assess the performance by varying the degrees of the local polynomials. The corresponding plots (See Figure 8) indicate that use of lower-order filters achieves a fit comparable to one by a relatively higher-degree polynomial. Since an actual waveform is not observed in practice and instead, we estimate it by optimizing over the observed labels via learning a graph filter, we theoretically show that the family of filters that we learn is a strict superset of the polynomial filter family. The same result holds for the families of the resulting adapted graphs. Theorem 4.2. Define H := {h(·) | ∀ possible K-degree polynomial parameterizations of h} to be the set of all K-degree polynomial filters, whose arguments are n × n diagonal matrices, such that a filter response over some Λ is given by h(Λ) for h(·) ∈ H. Similarly H′ := {h̃(·) | ∀ possible polynomial parameterizations of h̃} is set of all filters learn-able via PP-GNN , with h̃ = h + hf1 + hf2 , where h is parameterized by a K-degree polynomial supported over entire spectrum, hf1 and hf2 are adaptive filters parameterized by independent K
′-degree polynomials which only act on top and bottom t diagonal elements respectively, with t < n/2 and K ′ ≤ K; then H and H′ form a vector space, with H ⊂ H′. Also, dim(H
′) dim(H) = K+2K′+3 K+1 .
Corollary 4.2.1. The corresponding adapted graph families G := {UTh(·)UT | ∀h(·) ∈ H} and G′ := {UT h̃(·)UT | ∀h̃(·) ∈ H′} for any unitary matrix U form a vector space, with G ⊂ G′ and dim(G′) dim(G) = K+2K′+3 K+1 .
This implies that our model is learning from a more diverse space of filters and the corresponding adapted graphs. Moreover, the dimension of the space increases significantly by using just two adaptive filters. Note that learning from such a diverse region is feasible. This observation is possible from the proofs of Theorem 4.2 (A.4.4) and Corollary 4.2.1 (A.4.5). Using the adaptive filters without any filter with the entire spectrum as support results in learning a set of adapted graphs, Ĝ. This set is disjoint from G, with G′ = G ⊕ Ĝ. We conduct various ablative studies where we demonstrate the effectiveness of learning from Ĝ and G′.

heading: 5 EXPERIMENTS
text: We conduct comprehensive experiments to demonstrate the effectiveness and competitiveness of the proposed PP-GNN model over state-of-the-art (SOTA) models/methods and answer the following research questions.
1. [RQ1] Does the PP-GNN model outperform SOTA models/methods over a broad range of datasets (homophilic/heterophilic)?
2. [RQ2] Does the proposed model learn better and complex frequency responses that improve performance [Frequency response plots]?
3. [RQ3] With the PP-GNN model comprising of multiple filters, how does each filter contribute to achieving improved performance?
4. [RQ4 and RQ5] Does the model learn better embedding, and do we need a large number of eigenvalues/vectors and polynomials to get improved performance?
5. [RQ6] How does the training time of the PP-GNN and GPR-GNN models compare?
Datasets and Setup. We evaluate our model on several real-world heterophilic and homophilic datasets including a few large graphs for the node classification task. Detailed statistics of the benchmark datasets are provided in the appendix (A.5.1). The heterophilic datasets include Texas2, Cornell2, Wisconsin2, Chameleon, Squirrel, (Rozemberczki et al., 2021) and Flickr. We use the 10 random splits (48%/32%/20% of nodes for train/validation/test set) from Pei et al. (2020). For Flickr, Ogbn-Arxiv, and Wiki-CS, we use the data splits from Kim & Oh (2021). We also
2http://www.cs.cmu.edu/afs/cs.cmu.edu/project/theo-11/www/wwkb
evaluate on 8 homophilic datasets: Cora-Full, Ogbn-Arxiv, Wiki-CS, Citeseer, Pubmed, Cora, Computer, and Photos borrowed from Kim & Oh (2021). For the remaining homophilic datasets, we create 10 random splits for each dataset following Kim & Oh (2021). We use PCA node features (Kim & Oh (2021)) for the Chameleon and Cora-Full datasets. We report mean and standard deviation of test accuracy over splits to compare model performance.
Baselines. We compare our model against (a) filtering based approaches, GPR-GNN, FAGCN, APPNP, and LGC, (b) other state-of-the-art models include SUPERGAT, TDGNN, H2GCN, and GEOM-GCN, and (c) standard baselines, namely, LR (Logistic Regression), MLP (Multi Layer Perceptron), GCN, and SGCN. Detailed descriptions for all the baselines, hardware and software specifications are in the appendix (A.5.2 and A.5.5)
Hyperparameter Tuning. For the PP-GNN model, we separately partition the low-end and highend eigenvalues into several contiguous partitions and use shared filter parameters for frequencies of each partition. The number of partitions, which can be interpreted as the number of filters, is swept in the range [2,3,4,5,10,20]. The polynomial filter order is swept in the range [1,10] in steps of size 1. The number of eigenvalues/vectors are swept in the range [32, 64, 128, 256, 512, 1024]. In our experiments, we set ηl = ηh and we vary them in range (0, 1) and set ηgpr = 1 − ηl. We did hyperparameter tuning for other models as suggested in respective references and public repositories. More details are in the appendix (A.5.5).

heading: 5.1 RQ1: PP-GNN VERSUS SOTA MODELS
text: We present several important observations from Tables 1 and 2. We observe that the PP-GNN model consistently outperforms all models, including recent filtering approach based models on all datasets (both heterophilic and homophilic) except Pubmed and Cornell, where it achieves similar performance. This result demonstrates the effectiveness and robustness of our model across a wide variety of datasets. Furthermore, compared with the GPR-GNN model, learning piece-wise polynomial filters improves performance significantly over learning a single polynomial filter. In particular, PP-GNN achieves performance improvements of around 10% and 5% on the Chameleon and Squirrel datasets.

heading: 5.2 RQ2: ADAPTABLE FREQUENCY RESPONSES
text: We computed the frequency responses (i.e., h(λ)) of learned polynomials of the PP-GNN and GPRGNN models on several datasets, including Squirrel and Citeseer datasets shown in Figure 2. We observe from Figure 2a, though the GPR-GNN model can learn some variations at low/high frequencies, it is insufficient to achieve higher classification accuracy. On the other hand, the PPGNN model can capture complex shapes at the low and high ends of the spectrum, enabling it to achieve significantly improved test accuracy. We observed a similar phenomenon for the chameleon dataset as well. To illustrate another behaviour, we present the frequency responses for the Citeseer dataset in Figure 2b, and the responses are similar except for some variations at the low-end of the spectrum. Note that the GPR-GNN model does well on several homophilic datasets. These observations suggest that the PP-GNN model adapts very well in learning desired frequency responses, as dictated by the task at hand. We can observe such a behaviour for two other datasets in the appendix (A.5.6).

heading: 5.3 RQ3: PERFORMANCE COMPARISON OF DIFFERENT FILTERS
text: Recall that the PP-GNN model is a sum of polynomials model comprising of polynomial filters operating at different parts of the spectrum. Here, we study the importance and effect of using a combination of filters operating on different regions. For reference, we also compare the performance of GPR-GNN model that operates on the entire spectrum using a single filter. Several interesting observations are in order. From Table 3, we see that the best performance is achievable using high-frequency signals alone for the heterophilic datasets (Squirrel and Chameleon), suggesting that significant discriminatory information is available at high frequencies compared to low-frequency signals. In contrast, homophilic datasets (Cora and Citeseer) exhibit a reverse trend. On comparing the first (second) and third (fourth) row results, we see that having the GPR-GNN filter as part of the PP-GNN filter helps to get improved performance over individual filters (PP-GNN(Low) or PP-GNN (High)). Finally, the PP-GNN model can adapt well, capture contrasting information bands across datasets, and outperform the GPR-GNN model.

heading: 5.4 RQ4 - RQ6: VARYING NO. OF EVS, LEARNED EMBEDDING AND TIMING COMPARISON
text: We conducted ablative studies to see the effect of varying different hyperparameters of the PP-GNN model. From Figure 3a, we see that as small as 32 eigen components are sufficient to achieve nearbest performance on the homophilic datasets (Cora, Citeseer). On the other hand, peak performance is achieved for some intermediate (≈250-500) number of eigen components for the Squirrel and Chameleon datasets. Also, we note that we can achieve comparable performance to state of the art models on large datasets (Flickr and OGBN-Arxiv) using only 1024 components (See Tables 1 and 2). This study suggests that the PP-GNN model is scalable at least for medium sized graphs ( 100K nodes). More evaluation and investigation are needed to study performance versus computation cost trade-off for very large graphs (millions of nodes).
To assess the quality of learned embedding, we also created t-SNE plots and made a visual inspection. From Figure 3b and 3c, we see that the clusters belonging to different classes of the Squirrel dataset are better separated for the PP-GNN model compared to the GPR-GNN model. We also measured the time taken by the GPR-GNN and PP-GNN models. We observed that the GPRGNN model is only (1.2 − 2)× faster than the PP-GNN model (using 1024 eigencomponents) on many datasets. Thus, the extra training cost incurred by the PP-GNN model is not significantly different for practical purposes. More details are in the appendix (A.6).

heading: 6 CONCLUSION
text: Several proposed models attempt to be robust to the correlations between graph and node labels. We build on the filter-based approach of GPR-GNN (Chien et al., 2021). This work proposed an effective polynomial filter design. We combine GPR-GNN with additional polynomials that adapt specifically to low and high-frequency components. Our experiments demonstrate that such an approach can learn filter functions that improve performance on the node classification task. Our plots of these filter functions suggest that they are of high order on several datasets. It would be interesting to analyze these filter plots and identify some common characteristics. These will enable us to a] characterize the correlation of graphs and labels, b] further improve the performance and c] build robust graph privacy models. We plan to do this as our future work.

heading: 7 REPRODUCIBILITY STATEMENT
text: We have taken several steps to ensure the reproducibility of our work. First, we provide a detailed description of our model in Section 4. Section 5 provides a brief description of the datasets utilized in our experiments. A detailed description of the datasets, including their sources, statistics are given in the appendix (A.5.1). We also provide details on the number of splits and the method used to generate them. We also list out all the baselines utilized in the experiments in Section 5. A detailed description of these baselines and the parameter sweep ranges are given in the appendix (A.5.2). We also provide implementation detail of our model in the appendix (A.5.5). Additionally, we share an anonymized URL of the code for our model, and the data splits in the appendix (A.1).

heading: 8 ETHICS STATEMENT
text: Social network graphs form the most exciting application for GNNs. GNNs can be used to reveal user information that the user otherwise would have preferred to keep private. Privacy-preserving methods might rely on obfuscating the graph by adding spurious edges, and these edges will, in turn, reduce the graph’s correlation with the node labels. Earlier GNNs would likely not have predicted the target label with reasonable accuracy under this setting. Thereby, such obfuscation methods would give some level of privacy for the users. However, our model exhibits good performance even in the presence of a lower correlation. Such methods can make simple obfuscation based privacy approaches obsolete while still revealing important information about the users. On the other hand, our model can also give insights into what makes graphs reveal certain information and develop a more robust privacy model for graphs.

heading: A APPENDIX
text: The appendix is structured as follows. In Section A.1, we provide a URL with our code and data, along with some information about the contents in the URL. In Section A.2, we present additional evidence of the limitations of GPR-GNN. In Section A.3, we show a representative experiment that motivates Section 4. In Section A.4, we provide proofs for theorems and corollaries defined in Section 4.3. In Section A.5, we provide more details regarding the baselines, datasets and their respective splits. We also provide implementation details and rank PP-GNN against current SoTA. In Section A.6, we perform a timing analysis where we compare PP-GNN to GPR-GNN.

heading: A.1 CODE AND DATASET
text: Our code, along with the datasets and their respective splits, is available at https://tinyurl. com/PPGNN. A README file is available explaining how the code can be executed on various datasets. The requirements.txt file lists all the necessary python packages required to execute the code. The src folder also contains bash scripts to run the code on various datasets.

heading: A.2 MOTIVATION

heading: A.2.1 NODE FEATURE INDISTINGUISHABLY PLOTS
text: In the main paper (Figure 1a), we plot the average of pairwise distances between node features for the Cora dataset, after computing ÃjX for increasing j values, and showed that the mean pairwise node feature distance decreases as j increases. We observe that this is consistent across three more datasets: Citeseer, Chameleon and Squirrel. This is observed in Figure 4.
We also observed the mean of the variance of each dimension of node features, after computing ÃjX , for increasing j values. We observe that this mean does indeed reduce as the number of hops increase. We also observe that the variance of each dimension of node features reduces for Cora, Squirrel and Chameleon as the number of hops increase; however, we don’t observe such an explicit phenomenon for Citeseer. See Figure 5.

heading: A.2.2 EFFECT OF VARYING THE ORDER OF THE GPR-GNN POLYNOMIAL
text: In the main paper (Figure 1a), we plot the test accuracies of the GPR-GNN model while increasing the order of the polynomials for the Cora and Chameleon dataset, respectively. We observe that on increasing the polynomial order, the accuracies do not increase any further. We can show a similar phenomenon on two other datasets, Squirrel and Citeseer, in Figure 6.
In Section 3 of the main paper, we claim that due to the over-smoothing effect, even on increasing the order of the polynomial, there is no improvement in the test accuracy. Moreover, in Figure 2 we can see that our model can learn a complicated filter polynomial while GPR-GNN cannot. This section shows that even on increasing the order of the GPR-GNNpolynomial, neither does the test accuracy increase nor does the waveform become as complicated as PP-GNN. See Figure 7.

heading: A.3 FICTITIOUS POLYNOMIAL
text: In Section 4, we claim that having multiple disjoint low order polynomials can approximate a complicated waveform better than a higher-order polynomial. We create a representative experiment that shows this is indeed true by creating a fictitious complicated polynomial and trying to fit it using both a single polynomial and the disjoint multi polynomial. We observe in Figure 8 that the lower order disjoint polynomial fits the complicated waveform better than a higher-order single polynomial.

heading: A.4 PROPOSED APPROACH

heading: A.4.1 DETAILS REGARDING BOUNDARY REGULARIZATION
text: To induce smoothness in the learned filters, we add a regularization term that penalizes squared differences between function values of polynomials at knots (endpoints of contiguous bins). Our regularization term looks as follows:
m−1∑ i=1 exp−(σ max i −σ min i+1 ) 2 (hi(σ max i )− hi+1(σmini+1 ))2 (5)
In equation 5, σmaxi and σ min i refer to the maximum and minimum eigenvalues in σi (Refer to Section 4). This regularization term is added to the Cross-Entropy loss. We perform experiments with this model and report the performance in Table 4. We observe that we are able to reach similar performance even without the presence of this regularization term. Therefore, majority of the results reported in our Main paper are without this regularization term.

heading: A.4.2 NOTATION USED
text: Vectors are denoted by lower case bold Roman letters such as x, and all vectors are assumed to be column vectors. In the paper, h with any sub/super-script refers to a frequency response, which is
also considered to be a vector. A superscript T denotes the transpose of a matrix or vector; Matrices are denoted by bold Roman upper case letters, such as M. A field is represented by K; sets of real and complex numbers are denoted by R and C respectively. K[x1, . . . , xn] denotes a multivariate polynomial ring over the field K, in indeterminates x1, . . . , xn. Set of n × n square matrices with entries from some set S are denoted by Mn(S). Moore-Penrose pseudoinverse of a matrix A is denoted by A†. Eigenvalues of a matrix are denoted by λ, with λ1, λ2, . . . denoting a decreasing order when the eigenvalues are real. A matrix Λ denotes a diagonal matrix of eigenvalues. Set of all eigenvalues, i.e., spectrum, of a matrix is denoted by σA or simply σ when the context is clear. Lp norms are denoted by ‖ · ‖p. Frobenius norm over matrices is denoted by ‖ · ‖F . Norms without a subscript default to L2 norms for vector arguments and Frobenius norm for matrices. ⊕ denotes a direct sum. For maps fi defined from the vector spaces V1, · · · , Vm, with a map of the form f : V 7→W , with V = V1⊕V2⊕· · ·⊕Vm, the the phrase ”f : V 7→W by mapping f(vi) to fi(g(vi))” means that f maps a vector v = v1 + . . .+ vm with vi ∈ Vi to f1(g(v1)) + . . .+ fm(g(vm))

heading: A.4.3 PROOF OF THEOREM 4.1
text: Theorem. For any desired frequency response h∗, and an integer K ∈ N, let h̃ := h + hf , with hf having a continuous support over a subset of the spectrum, σf . Assuming h and hf to be parameterized by independent K and K ′-order polynomials p and pf respectively, with K ′ ≤ K, then there exists h̃, such that min ‖h̃ − h∗‖2 ≤ min ‖h − h∗‖2, where the minimum is taken over the polynomial parameterizations. Moreover, for multiple polynomial adaptive filters hf1 , hf2 , ..., hfm parameterized by independentK ′-degree polynomials withK ′ ≤ K but having disjoint, contiguous supports, the same inequality holds for h̃ = h+ ∑m i=1 hfi .
Proof. We make the following simplifying assumptions:
1. |σfi | > K, ∀i ∈ [m], i.e., that is all support sizes are lower bounded by K (and hence K ′)
2. All eigenvalues of the reference matrix are distinct
For methods that use a single polynomial filter, the polynomial graph filter, hK(Λ) = diag(Vγ) where γ is a vector of coefficients (i.e, γ parameterizes h), with eigenvalues sorted in descending order in components, and V is a Vandermonde matrix:
V =  1 λ1 λ 2 1 · · · λK1 1 λ2 λ 2 2 · · · λK2 ... ... ... . . .
... 1 λn λ 2 n · · · λKn  And to approximate a frequency response h∗, we have the following objective:
min ‖h− h∗‖22 := min γ ‖diag(h∗)− diag(Vγ)‖2F = minγ ‖h ∗ −Vγ‖22 = minγ ‖ep(γ)‖ 2 2
Where ‖‖F and ‖‖2 are the Frobenius and L2 norms respectively. Due to the assumptions, the system of equations h∗ = Vγ is well-defined and has a unique minimizer, γ∗ = V†h∗, and thus
‖ep(γ∗)‖ = minγ ‖ep(γ)‖. Next we break this error vector as:
ep(γ ∗) := h∗ −Vγ∗
= m∑ i=1 (h∗i −Viγ∗) + (h∗L −VLγ∗)
:= m∑ i=1 e∗pi + e ∗ pL
Where e∗pi := (h ∗ i −Viγ∗) with similar definition for epL ; h∗i is a vector whose value at components corresponding to the set σ(hfi) is same as that of h ∗ and rest are zero. Similarly, V∗i is a matrix whose rows corresponding to the set σ(hfi) are same as that of V with other rows being zero. Also, VL = V − ∑m i=0 Vi and h ∗ L = h ∗ − ∑m i=0 h ∗ i . Note that as a result of this construction, [ep∗i ] ∪ ep∗L is a linearly independent set since the supports [σ(hfi)] form a disjoint set (note the theorem statement). We split the proof in two cases:
Case 1: K ′ = K. We now analyze the case where we have m polynomial adaptive filters added, all having an order of K, where the objective is min ‖h̃− h∗‖, which can be written as:
min γ,[γi] ∥∥∥∥∥diag(h∗)− diag ( Vγ + m∑ i=0 Viγi )∥∥∥∥∥ 2
F
= min γ,[γi] ∥∥∥∥∥h∗ −Vγ − m∑ i=0 Viγi ∥∥∥∥∥ 2
2
= min γ,[γi]
‖eg(γ, [γi])‖22
Before characterizing the above system, we break a general error vector as:
eg(γ, [γi]) := h ∗ −Vγ − m∑ i=0 Viγi
= m∑ i=1 (h∗i −Vi(γ + γi)) + (h∗L −VLγ)
:= m∑ i=1 egi + egL
Where egi := (h ∗ i −Vi(γ + γi)) with similar definition for egL . Clearly, the systems of equations, egi = 0, ∀i and egL = 0 are well-defined due to the assumptions 1 and 2. Since all the systems of equations have independent argument, unlike in the polynomial filter case where the optimization is constrained over a single variable; one can now resort to individual minimization of squared norms of egi which results in a minimum squared norm of eg . Thus, we can set:
γ = V†Lh ∗ L = γ ∗ g γi = V † ih ∗ i −V † Lh ∗ L = γ ∗ i , ∀i ∈ [m]
to minimize squared norms of egi and egL . Note that [egi ] ∪ egL is a linearly independent set since the supports [σ(hfi)] form a disjoint set and by the above construction, this is also an orthogonal set, and hence we have ‖eg‖2 = ∑m i=1 ‖egi‖ 2 + ‖egL‖ 2, and hence the above assignment implies:
‖eg(γ∗g , [γ∗i ])‖ = min γ,[γi] ‖eg(γg, [γi])‖ := min ‖h̃− h∗‖2
Hence, it follows that, minx ‖h∗i −Vix‖2 = ∥∥e∗gi∥∥2 ≤ ∥∥e∗pi∥∥2 = ‖h∗i −Viγ∗‖2 and minx ‖h∗L −
VLx‖2 = ∥∥e∗gL∥∥2 ≤ ∥∥e∗pL∥∥2 = ‖h∗L −VLγ∗‖2. Hence,
m∑ i=1 ∥∥e∗gi∥∥2 + ∥∥e∗gL∥∥2 ≤ m∑ i=1 ∥∥e∗pi∥∥2 + ∥∥e∗pL∥∥2 min ‖h̃− h∗‖ ≤ min ‖h− h∗‖
Case 2: K ′ < K. We demonstrate the inequality showing the existence of an h̃ that achieves a better approximation error. By definition, the minimum error too will be bounded above by this
error. For this, we fix γ, the parameterization of h as γ = V†h∗ = γ∗p (say). Note that γp∗ = arg minγ ‖ep(γ)‖. Now our objective function becomes
eg(γ ∗ p , [γi]) := h ∗ −Vγ∗p − m∑ i=0 V′iγi
= m∑ i=1 (h∗i −Viγ∗p + V′iγi) + (h∗L −VLγ∗p)
= m∑ i=1 e′gi + e ′ gL
Where h∗i , h ∗ L,Vi,VL have same definitions as that in case 1 and V ′ i is a matrix containing first K ′ + 1 columns of Vi as its columns (and hence has full column rank), and, γi ∈ RK ′+1. By this construction, we have ‖eg(γ∗p ,0)‖ = min
γ ‖ep(γ)‖ = ‖ep(γ∗p)‖
Our optimization objective becomes min[γi] ‖eg(γ∗p , [γi])‖, which is easy since the problem is wellposed by assumption 1 and 2. The unique minimizer of this is obtained by setting
γi = V ′† i (h ∗ i −Viγ∗p) = γ∗i (say) ∀i ∈ [m]
Now, ‖eg(γ∗p , [γ∗i ])‖ = min
[γi] ‖eg(γ∗p , [γi])‖ ≤ ‖eg(γ∗p ,0)‖ = min γ ‖ep(γ)‖ = min ‖h− h∗‖

heading: By the definition of minima, minγ,[γi] ‖eg(γ, [γi])‖ ≤ min[γi] ‖eg(γ∗p , [γi])‖, and by the definition,
text: min ‖h̃− h∗‖ = minγ,[γi] ‖eg(γ, [γi])‖, we have:
min ‖h̃− h∗‖ ≤ min ‖h− h∗‖

heading: A.4.4 PROOF OF THEOREM 4.2
text: Theorem. Define H := {h(·) | ∀ possible K-degree polynomial parameterizations of h} to be the set of all K-degree polynomial filters, whose arguments are n × n diagonal matrices, such that a filter response over some Λ is given by h(Λ) for h(·) ∈ H. Similarly H′ := {h̃(·) | ∀ possible polynomial parameterizations of h̃} is set of all filters learn-able via PP-GNN , with h̃ = h + hf1 + hf2 , where h is parameterized by a K-degree polynomial supported over entire spectrum, hf1 and hf2 are adaptive filters parameterized by independent K
′-degree polynomials which only act on top and bottom t diagonal elements respectively, with t < n/2 and K ′ ≤ K; then H and H′ form a vector space, with H ⊂ H′. Also, dim(H
′) dim(H) = K+2K′+3 K+1 .
Proof. We start by constructing the abstract spaces on top of the polynomial vector space. Consider the set of all the univariate polynomials having degree at most K in the vector space over the ring Kxn := K[x1, . . . , xn] where K is the field of real numbers. Partition this set into n subsets, say V1, . . . , Vn, such that for i ∈ [n], Vi contains all polynomials of degree up to K in xi. It is easy to see that V1, . . . , Vn are subspaces of K[x1, . . . , xn]. Define V = V1 ⊕ V2 ⊕ · · · ⊕ Vn where ⊕ denotes a direct sum. Define the matrix Di[c] whose (i, i)th entry is c and all the other entries are zero. For i ∈ [n], define linear maps φi : Vi → Mn (Kxn) by f(xi) 7→ Di[f(xi)]. Im(φi) forms a vector space of all diagonal matrices, whose (i, i) entry is the an element of Vi. Generate a linear map φ : V → Mn(Kxn) by mapping φ(f(xi)) to φi(f(xi)) for all i ∈ [n] as the components of the direct sum present in its argument. Note that φi for i ∈ [n] are injective maps, making φ an injective map. This implies that H ⊂ Im(φ) is a subspace with basis Bh := {φ(x01 + · · ·+ x0n), φ(x1 + · · ·+ xn), . . . , φ(x K 1 + · · ·+xKn )}, making dim(H) = K+1. Similarly we have, H′ ⊂ Im(φ), a subspace
with basis Bh′ := Bh ⋃ {φ(x01 + · · ·+x0t +0+ · · ·+0), φ(x1 + · · ·+xt+0+ · · ·+0), . . . , φ(xK ′ 1 +
· · ·+xK′t + 0 + · · ·+ 0)} ⋃ {φ(0 + · · ·+ 0 +x0n−t+1 + · · ·+x0n), φ(0 + · · ·+ 0 +xn−t+1 + · · ·+ xn), . . . , φ(0+ · · ·+0+xK ′ n−t+1 + · · ·+xK ′
n )} where x0i and 0 are the corresponding multiplicative and additive identities of Kxn, implying H ⊂ H′ and dim(H′) = K + 2K ′ + 1.

heading: A.4.5 PROOF OF COROLLARY 4.2.1
text: Corollary. The corresponding adapted graph families G := {UTh(·)UT | ∀h(·) ∈ H} and G′ := {UT h̃(·)UT | ∀h̃(·) ∈ H′} for any unitary matrix U form a vector space, with G ⊂ G′ and dim(G′) dim(G) = K+2K′+3 K+1 .
Proof. Consider the injective linear maps f1, f2 : Mn(Kxn) → Mn(Kxn) as f1(A) = UTA and f2(A) = AU. Define f3 : H → Mn(Kxn) and f4 : H′ → Mn(Kxn) as f3(A) = (f1 ◦ f2)(A) for A ∈ H and f4(A) = (f1 ◦ f2)(A) for A ∈ H′. Since U is given to be a unitary matrix, f3 and f4 are monomorphisms. Using this with the result from Theorem 4.2, H ⊂ H′, we have G ⊂ G′.

heading: A.5 EXPERIMENTS

heading: A.5.1 DATASETS
text: We evaluate on multiple benchmark datasets to show the effectiveness of our approach. Detailed statistics of the datasets used are provided in Table 5. We borrowed Texas, Cornell, Wisconsin from WebKB3, where nodes represent web pages and edges denote hyperlinks between them. Actor is a co-occurrence network borrowed from Tang et al. (2009), where nodes correspond to an actor, and and edge represents the co-occurrence on the same Wikipedia page. Chameleon, Squirrel are borrowed from Rozemberczki et al. (2021). Nodes correspond to web pages and edges capture mutual links between pages. For all benchmark datasets, we use feature vectors, class labels from Kim & Oh (2021). For datasets in (Texas, Wisconsin, Cornell, Chameleon, Squirrel, Actor), we use 10 random splits (48%/32%/20% of nodes for train/validation/test set) from Pei et al. (2020). We borrowed Cora, Citeseer, and Pubmed datasets and the corresponding train/val/test set splits from Pei et al. (2020). The remaining datasets were borrowed from Kim & Oh (2021). We follow the same dataset setup mentioned in Kim & Oh (2021) to create 10 random splits for each of these datasets. We also experiment with two slightly larger datasets Flickr Chua et al. (July 8-10, 2009) and OGBN-arXiv Hu et al. (2020). We use the publicly available splits for these datasets.

heading: A.5.2 BASELINES
text: We provide the methods in comparison along with the hyper-parameters ranges for each model. For all the models, we sweep the common hyper-parameters in the same ranges. Learning rate is swept over [0.001, 0.003, 0.005, 0.008, 0.01], dropout over [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8], weight decay over [1e− 4, 5e− 4, 1e− 3, 5e− 3, 1e− 2, 5e− 2, 1e− 1], and hidden dimensions over [16, 32, 64]. For model-specific hyper-parameters, we tune over author prescribed ranges. We use undirected graphs with symmetric normalization for all graph networks in comparison. For all models, we report test accuracy for the configuration that achieves the highest validation accuracy. We report standard deviation wherever applicable.
LR and MLP: We trained Logistic Regression classifier and Multi Layer Perceptron on the given node features. For MLP, we limit the number of hidden layers to one.
GCN: We use the GCN implementation provided by the authors of Chien et al. (2021). Note: We observe discrepancies in performance with various implementations of GCN. We make a note of this in Section A.7.
3http://www.cs.cmu.edu/afs/cs.cmu.edu/project/theo-11/www/wwkb
SGCN: SGCN (Wu et al., 2019) is a spectral method that models a low pass filter and uses a linear classifier. The number of layers in SGCN is treated as a hyper-parameter and swept over [1, 2].
SUPERGAT: SUPERGAT (Kim & Oh, 2021) is an improved graph attention model designed to also work with noisy graphs. SUPERGAT employs a link-prediction based self-supervised task to learn attention on edges. As suggested by the authors, on datasets with homophily levels lower than 0.2 we use SUPERGATSD. For other datasets, we use SUPERGATMX. We rely on authors code4 for our experiments.
GEOM-GCN: GEOM-GCN (Pei et al., 2020) proposes a geometric aggregation scheme that can capture structural information of nodes in neighborhoods and also capture long range dependencies. We quote author reported numbers for Geom-GCN. We could not run Geom-GCN on other benchmark datasets because of the unavailability of a pre-processing function that is not publicly available.
H2GCN: H2GCN (Zhu et al., 2020) proposes an architecture, specially for heterophilic settings, that incorporates three design choices: i) ego and neighbor-embedding separation, higher-order neighborhoods, and combining intermediate representations. We quote author reported numbers where available, and sweep over author prescribed hyper-parameters for reporting results on the rest datasets. We rely on author’s code5 for our experiments.
FAGCN: FAGCN (Bo et al., 2021) adaptively aggregates different low-frequency and highfrequency signals from neighbors belonging to same and different classes to learn better node representations. We rely on author’s code6 for our experiments.
APPNP: APPNP (Klicpera et al., 2019) is an improved message propagation scheme derived from personalized PageRank. APPNP’s addition of probability of teleporting back to root node permits it to use more propagation steps without oversmoothing. We use GPR-GNN’s implementation of APPNP for our experiments.
LIGHTGCN: LIGHTGCN (Navarin et al., 2020) is a spectrally grounded GCN that adapts the entire eigen spectrum of the graph to obtain better node feature representations.
GPR-GNN: GPR-GNN (Chien et al., 2021) adaptively learns weights to jointly optimize node representations and the level of information to be extracted from graph topology. We rely on author’s code7 for our experiments.
TDGNN: TDGNN (Wang & Derr, 2021) is a tree decomposition method which mitigates feature smoothening and disentangles neighbourhoods in different layers. We rely on author’s code8 for our experiments.
ARMA: ARMA (Bianchi et al., 2021) is a spectral method that uses K stacks of ARMA1 filters in order to create an ARMAK filter (an ARMA filter of order K). Since (Bianchi et al., 2021) do not specify a hyperparameter range in their work, following are the ranges we have followed: GCS stacks (S): [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], stacks’ depth(T ): [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. However we only select configurations such that the number of learnable parameters are less than or equal to those in PP-GNN. The input to the ARMAConv layer are the node features and the output is the number of classes. This output is then passed through a softmax layer. We use the implementation from the official PyTorch Geometric Library 9
BernNet: BernNet (He et al., 2021) is a method that approximates any filter over the normalised Laplacian spectrum of a graph, by a Kth Order Bernstein Polynomial Approximation. We use the model specific hyper-parameters prescribed by the authors of the paper. We vary the Propagation Layer Learning Rate as follows: [0.001, 0.002, 0.01, 0.05]. We also vary the Propagation Layer
4https://github.com/dongkwan-kim/SuperGAT 5https://github.com/GemsLab/H2GCN 6https://github.com/bdy9527/FAGCN 7https://github.com/jianhao2016/GPRGNN 8https://github.com/YuWVandy/TDGNN 9https://pytorch-geometric.readthedocs.io/en/latest/_modules/torch_
geometric/nn/conv/arma_conv.html#ARMAConv
Dropout as follows: [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8]. We rely on the authors code 10 for our experiments.
AdaGNN: AdaGNN (Dong et al., 2021) is a method that captures the different importance’s for varying frequency components for node representation learning. We use the model specific hyperparameters prescribed by the authors of the paper. The No. of Layers hyper-parameter is varied as follows: [2, 4, 8, 16, 32, 128]. We rely on the authors code11 for our experiments.

heading: A.5.3 COMPARISON AGAINST ADDITIONAL BASELINES
text: We performed experiments over the same datasets using three new baselines: ARMA (Bianchi et al. (2021)), BernNet (He et al. (2021)) and AdaGNN (Dong et al. (2021)). Details about these baselines are given in Section A.5.2. The results are reported in Tables 6 and 7. The results are obtained after thorough optimization over the respective hyperparameters, the details of which are listed in A.5.2:
For AdaGNN and BernNet, we use the authors code and tune it over the hyperparameters as provided in the paper. For ARMA, we use the official PyTorch Geometric implementation (see A.5.2). As a sanity check, we also tested ARMA on the node classification datasets described in the paper and were able to reproduce similar numbers.

heading: A.5.4 COMPARISON AGAINST GENERAL FIR FILTERS
text: Instead of using a polynomial filter, we can use a general FIR filter (GFIR) which is described by the following equation:
Z = K∑ k=0 SkXHk
where S is the graph shift operator (which in our case is Ã), X is the node feature matrix and Hk’s are learnable filter matrices. One can see GCN, SGCN, GPR-GNN as special cases of this GFIR filter, which constrain the Hk in different ways.
We first demonstrate that constraint on the GFIR filter is necessary for getting improvement in performance, particularly on heterophilic datasets. Towards this, we build two versions of GFIR: one with regularization (constrained), and the other without regularization (unconstrained). We ensure that that the number of trainable parameters in these models are comparable to those used in PP-GNN. We provide further details of the versions of the GFIR models below and report the results in Table 8 below:
10https://github.com/ivam-he/BernNet 11https://github.com/yushundong/AdaGNN
• Unconstrained Setting: In this setting, we do not impose any regularization constraints such as dropout and L2 regularization.
• Constrained Setting: In this setting, we impose dropouts as well as L2 regularization on the GFIR model. Both dropouts and L2 regularization were applied on the Hk’s (the learnable filter matrices from the above equation).
We also compare PP-GNN (the proposed model) as well as GPR-GNN to the General FIR filter model (GFIR).

heading: Train Acc /Test Acc Computers Chameleon Citeseer Cora Squirrel Texas Wisconsin
text: We can make the following observation from the results reported in Table 8:
• Firstly, constrained GFIR performs better than the unconstrained version, with performance lifts of up to ∼10%. This suggests that regularization is important for GFIR models.
• GPR-GNN outperforms the constrained GFIR version. It is to be noted that GPR-GNN further restricts the space of graphs explored as compared to GFIR. This suggests that regularization beyond simple L2/dropout kind of regularization (polynomial filter) is beneficial.
• PP-GNN performs better than GPR-GNN. Our model slightly expands the space of graphs explored (as compared to GPR-GNN, but lesser than GFIR), while retaining good performance. This suggests that there is still room for improvement on how regularization is done.
PP-GNN has shown one possible way to constrain the space of graphs while improving performance on several datasets, however, it remains to be seen whether there are alternative methods that can do even better. We hope to study and analyze this aspect in the future.
A.5.5 IMPLEMENTATION DETAILS
In this subsection, we present several important points that are useful for practical implementation of our proposed method and other experiments related details. Our approach is based on adaptation of eigen graphs constructed using eigen components. Following Kipf & Welling (2017), we use a symmetric normalized version (Ã) of adjacency matrix A with self-loops: Ã = D̃− 1 2 (A + I)D̃− 1 2
where D̃ii = 1 + Dii, Dii = ∑ j Aij and D̃ij = 0, i 6= j. We work with eigen matrix and eigen
values of Ã.
To reduce the learnable hyper-parameters, we separately partition the low-end and high-end eigen values into several contiguous bins and use shared filter parameters for each of these bins. The number of bins, which can be interpreted as number of filters, is swept in the range [2, 3, 4, 5, 10, 20]. The orders of the polynomial filters are swept in the range [1, 10] in steps of 1. The number of EVD components are swept in the range [32, 64, 128, 256, 512, 1024]. In our experiments, we set ηl = ηh and we vary the ηl parameter in range (0, 1) and ηgpr = 1− ηl. For optimization, we use the Adam optimizer (Kingma & Ba, 2015). We set early stopping to 200 and the maximum number of epochs to 1000. We utilize learning rate with decay, with decay factor set to 0.99 and decay frequency set to 50. All our experiments were performed on a machine with Intel Xeon 2.60GHz processor, 112GB Ram, Nvidia Tesla P-100 GPU with 16GB of memory, Python 3.6, and PyTorch 1.9.0 (Paszke et al., 2019). We used Optuna (Akiba et al., 2019) to optimize the hyperparameter search.

heading: A.5.6 ADAPTABLE FREQUENCY RESPONSES
text: In Figure 2 of the main paper, we observe that PP-GNN learns a complicated frequency response for a heterophilic dataset (Squirrel) and a simpler frequency response for a homophilic dataset (Citeseer). We observe that this trend follows for two other datasets Chameleon (heterophilic) and Computer (homophilic). See Figure 9.

heading: A.5.7 DOES THE POLYNOMIAL INITIALIZATION MATTER FOR GIVING GOOD PERFORMANCE?
text: In our experiments, we have varied the polynomial initialization schemes namely: [PPR, NPPR and Random], as described in GPRGNN paper Chien et al. (2021). To understand whether this polynomial initialization is important, we perform an ablative study to compare the various initialization schemes. We observe that our model (PP-GNN) retains similar performance across all initialization schemes. For brevity, we report results on six datasets (Cora, Citeseer, Squirrel, Chameleon, Texas and Wisconsin). The results can be found in Table 9. Our experiment results suggest that tuning polynomial initialization for each partition can indeed be avoided without compromising on the performance.

heading: A.5.8 DOES THE MLP EVEN MATTER?
text: In PP-GNN there is a two layered MLP (that transforms the input node features) followed by a single graph filtering layer similar to GPR-GNN.
To understand MLP’s significance, we ran an additional experiment, where we have used a single linear layer, instead of the two layered MLP. We continue to observe competitive (with respect to our original PP-GNN model) performance, across most datasets. The results can be found in Table 10. Also, the two layer MLP is not the significant contributor towards performance. This can also be seen by comparing GPR-GNN’s performance with that of LGC’s. LGC can be interpreted as a linear version of GPR-GNN, and achieves comparable performance as GPR-GNN.
Note that PP-GNN (Linear) was not thoroughly trained as the PP-GNN (Original) due to time constraints. The hyperparameter ranges are explained right below the next paragraph.
Table 10 above seems to suggest that PP-GNN (Linear) is competitive to PP-GNN (Original). We believe the difference in performance is an artifact of inadequate tuning. We can infer that adding MLP may give marginal improvements over the linear version. This phenomenon is also observed in GPR-GNN. To illustrate this, we can compare GPR-GNN with LGC (linear version of GPR-GNN). We can observe in the above table that GPR-GNN and LGC are comparable in performance.
For this additional set of experiments (Linear PP-GNN), due to time constraints, we have reduced the ranges for several hyper-parameters of our model, while also reducing the number of Optuna Trials to 50. These reduced ranges are as follows:
• Learning Rate: [0.003, 0.005] • Weight Decay : [0.0001, 0.001] • Dropout: [0.3, 0.5] • Hidden Dims: [32, 64] • Order of the Polynomial: [2, 4] • No. of Partitions (Buckets): [2, 4] • Number of Eigen Value: [256, 1024] • η: Sampled uniformly between [0, 1] using Optuna
Despite this smaller hyper-parameters search space, we are only shy by a few percentage points on the Texas dataset. We believe this difference can be recovered with thorough tuning.

heading: A.6 TRAINING TIME ANALYSIS
text: In the following subsections, we provide comprehensive timing analysis.
Computational Complexity: Listed below is the computational complexity for each piece in our model for a single forward pass. Notation n: number of nodes, |E|: the number of edges, A: symmetric normalized adjacency matrix, F : features dimensions, d: hidden layer dimension, C: number of classes, e∗ denotes the cost of EVD, K: polynomial order/hop order, l: number of eigenvalues/vectors in a single partition of spectrum (for implementation, we keep l same for all such intervals), m: number of partitions of a spectrum.
• MLP: O(nFd+ ndC) • GPR-term: O(K|E|C) + O(nKC). The first term is the cost for computing AKf(X) for
sparse A. The second term is the cost of summation ∑ k A
kf(X). • Excess terms for PP-GNN: O(mnlC). This is obtained by the optimal matrix multiplica-
tion present in Equation 3 of the main paper (Ui is n× l, Hi(γi) is l × l, Z0() is n× C). The additional factorm is because we havem different contiguous intervals/different polynomials. Note that typically n is much larger than l.
• EVD-term: e∗, the complexity for obtaining the eigenvalues/vectors of the adjacency matrix, which is usually very sparse for the observed graphs. Most publicly available solvers for this task utilize Lanczos’ algorithm (which is a specific case of a more general Arnoldi iteration). However, the convergence bound of this iterative procedure depends upon the starting vectors and the underlying spectrum (particularly the ratio of the absolute difference of two largest eigenvalues to the diameter of the spectrum) [Saad (1980), LI
(2010), Cullum & Willoughby (2002)]. Lanczos’ algorithm is shown to be a practically efficient way for obtaining extreme eigenpairs for a similar and even very large systems. We use ARPACK’s built-in implementation to precompute the eigenvalues/vectors for all datasets before training, thus amortizing this cost across training with different hyper-parameters configuration.
Per Component Timing Breakup: In Table 11, we provide a breakdown of cost incurred in seconds for different components of our model. Since the eigenpairs’ computation is a one time cost, we amortize this cost over the total hyper-parameters configurations and report the effective training time in the last column on of Table 11.
Average Training Time: In Table 12, we report the training time averaged over 100 hyperparameter configurations for several models. To understand the relative performance of our model with respect to GCN, we compute the relative time taken and report it in Table 13. We can observe in Table 13 that PP-GNN is ∼ 4x slower than GCN, ∼ 2X slower than GPR-GNN and BernNET, and ∼ 2X faster than AdaGNN. However, it is important to note that in our average training time, the time taken to compute K top and bottom eigenvalues/vectors is amortized across the number of trials
Can we cut down the total training cost for PP-GNN? To study this, we analyzed the chosen hyper-parameters for our 100 configurations experiments and prescribe below the new ranges for hyperparameters. Reduced hyperparameters space:
• Learning Rate: [0.003, 0.005]
• Weight Decay : [0.0001, 0.001]
• Dropout: [0.3, 0.5]
• Hidden Dims: [32, 64]
• Order of the Polynomial: [2, 4]
• No. of Partitions (Buckets): [2, 4]
• Number of Eigen Value: [256, 1024]
• η: Sampled uniformly between [0, 1] using Optuna
We further restrict the total number of Optuna trials to 20. Note that as shown in Section A.5.7, our model is insensitive to the polynomial initialization scheme. Hence we use ‘random’ initialization, thereby reducing the number of hyperparameters even further. The total time taken by PP-GNN (including the one-time cost for obtaining eigenpairs) for training (and performing hyperparameter optimization) over 20 hyper-parameter configurations is reported in Table 14. For comparison, we also show the total training time for the 100 configurations experiments in Table 15 which too includes the time taken to obtain min(n, 1024) eigenpairs. From Tables [14,15] we can observe that the total training time has reduced by ∼ 5x. In Table 16, we can observe that even with the reduced hyper-parameter range, our model performs competitively or better.

heading: A.7 DISCREPANCIES IN PERFORMANCE BETWEEN DIFFERENT GCN IMPLEMENTATIONS
text: We tested three different implementations of GCN(s) (Kipf & Welling, 2017) provided by 1) Thomas Kipf and Max Welling 12 (In TensorFlow), 2) The implementation13 of the authors of the GPR-GNN paper (Chien et al., 2021) (in PyTorch), 3) The code base 14 publicly available at the OGBN-ArXiv leader-board 15 (in PyTorch), on the OGBN-ArXiv dataset.
12https://github.com/tkipf/gcn 13https://github.com/jianhao2016/GPRGNN 14https://github.com/snap-stanford/ogb/tree/master/examples/
nodeproppred/arxiv 15https://ogb.stanford.edu/docs/leader_nodeprop/#ogbn-arxiv
Note that the authors of the leader-board code report a test accuracy of 71.74% on the OGBN-ArXiV dataset. Extensive tuning was done on the Kipf and Welling implementation in order to match the results of the leaderboard code. However, even after extensive hyper-parameter tuning, the test accuracy increased to only 65.53%, which is still much lesser.
Two other implementations that we tested were: a) The GCN implementation from the authors of the GPR-GNN paper (Chien et al., 2021) (in PyTorch) b) The leader-board implementation (in PyTorch). The results for all the three implementations on the OGBN-ArXiv Dataset can be found in Table 17.

heading: Dataset Chameleon Citeseer Computers Cora Cora-Full Photo Pubmed Squirrel Texas Wisconsin OGBN-ArXiv
text: Besides different implementations of GCN(s), we find that batch normalization is done in the implementation used in the GCN leader-board code. We see that (in Table 17) the performance numbers vary depending on the implementation. Also, batch normalization helps to get nearly 2% improvement for this dataset. Since we did not perform batch normalization in all our experiments we find a performance gap of ∼ 2% for this dataset in results reported in our paper (against the leader-board code).
We have also compared the GCN implementation from Kipf and Welling as well as the GCN implementation from the authors of the GPR-GNN across various datasets in Table 18.
We can see that both the TensorFlow and the PyTorch implementations lead to a similar performance for most of the datasets. However, there are some differences in performance for some datasets (Wisconsin, Cora-Full (PCA), OGBN-ArXiv, Computer and Photos). We have ensured that hyperparameters, preprocessing steps etc., are same across both these implementations. We believe that these differences can be attributed to the different internal workings of TensorFlow and PyTorch.
A further study is required to understand where these gaps are coming from.
Note: In the main paper (in Tables [1, 2]), we report the results obtained on the GCN implementation by the authors of the GPR-GNN paper.

heading: Dataset Chameleon Citeseer Computers Cora Cora-Full Cornell Flickr Photo Pubmed Squirrel Texas Wikics Wisconsin OGBN-ArXiv

Review summaries:
This paper considers the problem of designing polynomial graph filters for use in graph neural networks. The motivation for doing so stems from the notions of homophily and heterophily in labeling the nodes of a graph: the former requiring the use of low-pass filters, and the latter requiring the use of high-pass filters. Here, the authors consider the design of spline polynomials for spectral-domain filtering, as opposed to the typical ``global'' polynomials used in graph filter design. That is, they first partition the spectrum of the graph matrix into low-pass, band-pass, and high-pass components, learn polynomials over each partition, and then combine these polynomials, enforcing continuity at the boundary of each interval via a penalty function.  Of course, doing this requires a full eigendecomposition of the graph matrix, which has high complexity. To ameliorate this issue, the authors propose an efficient variant of their spline spectral filter. This approach only uses a coarse partition of the spectrum, rather than using the full eigendecomposition.

This paper focuses on the task of node classification (semi-supervised learning) in heterophilic graphs. The authors identify that low-pass filter-based GNNs perform poorly on heterophilic graphs, precisely because they tend to smooth out differences in neighboring node features. Thus, they propose to look into filters capable of learning high-frequency content that can thus learn labels even if neighboring nodes do not share a label. In particular, the authors propose to learn different low-order polynomials for different parts of the spectrum.

 This paper proposes PP-GNN, a novel graph neural network that learns multiple adaptive polynomial filters acting on different subsets of the eigenvalues. The authors combine GPR-GNN with existing efficient algorithms for generating top and bottom eigen components to reduce the expensive complexity of eigendecomposition. They show that the piece-wise polynomial method can approximate a latent optimal filter better than a single polynomial in theory. 

In this paper, the authors aim to develop GNN that can better adapt to the given prediction task (both homophily and heterophiliy). Specifically, the authors extend the existing polynomial filter and propose to learn a filter function as a sum of polynomials over different subsets of the eigenvalues. The effectiveness of the proposed GNN architecture is demonstrated on diverse node classification tasks. The ablation studies were carried out to understand the proposed GNN architecture.