{'fraction_non_alphanumeric': 0.03725366098513782, 'fraction_numerical': 0.013552737827633618, 'mean_word_length': 4.705640642510272, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
NAFS: A Simple yet Tough-to-beat Baseline for Graph Representation Learning

Recently, graph neural networks (GNNs) have shown prominent performance in graph representation learning by leveraging knowledge from both graph structure and node features. However, most of them have two major limitations. First, GNNs can learn higher-order structural information by stacking more layers but can not deal with large depth due to the over-smoothing issue. Second, it is not easy to apply these methods on large graphs due to the expensive computation cost and high memory usage. In this paper, we present node-adaptive feature smoothing (NAFS), a simple non-parametric method that constructs node representations without parameter learning. NAFS first extracts the features of each node with its neighbors of different hops by feature smoothing, and then adaptively combines the smoothed features. Besides, the constructed node representation can further be enhanced by the ensemble of smoothed features extracted via different smoothing strategies. We conduct experiments on four benchmark datasets on two different application scenarios: node clustering and link prediction. Remarkably, NAFS with feature ensemble outperforms the state-of-the-art GNNs on these tasks and mitigates the aforementioned two limitations of most learning-based GNN counterparts.


Introduction
In recent years, graph representation learning has been extensively applied in various application scenarios, such as node clustering, link prediction, node classification, and graph classification (Kipf & Welling, 2016b;Wu et al., 2020a;1 Zhang et al., 2020;Miao et al., 2021a;Jiang et al., 2022;Wang et al., 2016;Wu et al., 2020b;Miao et al., 2021b). The goal of graph representation learning is to encode graph information to node embeddings. Traditional graph representation learning methods, such as DeepWalk (Perozzi et al., 2014), Node2vec (Grover & Leskovec, 2016), LINE (Tang et al., 2015), and ComE (Cavallari et al., 2017) merely focus on preserving graph structure information. GNNbased graph representation learning has attracted intensive interest by combining knowledge from both graph structure and node features. While most of these GNN-based methods are designed based on Graph AutoEncoder (GAE) and Variational Graph AutoEncoder (VGAE) (Kipf & Welling, 2016b), these methods share two major limitations: Shallow Architecture. Previous work shows that although stacking multiple GNN layers in Graph Convolutional Network (GCN) (Kipf & Welling, 2016a) is capable of exploiting deep structural information, applying a large number of GNN layers might lead to indistinguishable node embeddings, i.e., the over-smoothing issue (Li et al., 2018;Zhang et al., 2021a). Therefore, most state-of-the-art GNNs resort to shallow architectures, which hinders the model from capturing long-range dependencies.
Low Scalability. GNN-based graph representation learning methods can not scale well to large graphs due to the expensive computation cost and high memory usage. Most existing GNNs need to repeatedly perform the computationally expensive and recursive feature smoothing, which involves the participation of the entire graph at each training epoch. (Zhang et al., 2022) Furthermore, most methods adopt the same training loss function as GAE, which introduces high memory usage by storing the dense-form adjacency matrix on GPU. For a graph of size 200 million, its dense-form adjacency matrix requires a space of roughly 150GB, exceeding the memory capacity of the current powerful GPU devices.
To tackle these issues, we propose a new graph representation learning method, which is embarrassingly simple: just smooth the node features and then combine the smoothed features in a node-adaptive manner. We name this method node-adaptive feature smoothing (NAFS), and its goal is to construct better node embeddings that integrate the in-arXiv:2206.08583v1 [cs.LG] 17 Jun 2022 formation from both graph structural information and node features. Based on the observation that different nodes have highly diverse "smoothing speed", NAFS adaptively smooths each node feature and takes advantage of both low-order and high-order neighborhood information of each node. In addition, feature ensemble is also employed to combine the smoothed features extracted via different smoothing operators. Since NAFS is training-free, it significantly reduces the training cost and scales better to large graphs than most GNN-based graph representation learning methods. This paper is not meant to diminish the current advancements in GNN-based graph representation learning approaches. Instead, we aim to introduce an easier way to obtain high-quality node embeddings and understand the source of performance gains of these approaches better. Feature smoothing could be a promising direction towards a more simple and effective integration of information from both graph structure and node features.
Our contributions are as follows: (1) New perspective. To the best of our knowledge, we are the first to explore the possibility that simple feature smoothing without any trainable parameters could even outperform state-of-theart GNNs; this incredible finding opens up a new direction towards efficient and scalable graph representation learning.
(2) Novel method. We propose NAFS, a nodeadaptive feature smoothing approach along with various feature ensemble strategies, to fully exploit knowledge from both the graph structure and node features. (3) State-of-the-art performance. We evaluate the effectiveness and efficiency of NAFS on different datasets and graphbased tasks, including node clustering and link prediction. Empirical results demonstrate that NAFS performs comparably with or even outperforms the state-of-the-art GNNs, and achieves up to two orders of magnitude speedup. In particular, on PubMed, NAFS outperforms GAE (Kipf & Welling, 2016b) and AGE (Cui et al., 2020) by a margin of 9.0% and 3.8% in terms of NMI in node clustering, while achieving up to 65.4× and 88.6× training speedups, respectively.

Preliminary
In this section, we first explain the notations and problem formulation. Then, we review current GNNs and GNNbased graph representation learning.

Notations and Problem Formulation.
In this paper, we consider an undirected graph G = (V,E) with |V| = n nodes and |E| = m edges. Here we suppose that m ∝ n as it is the case in most real-world graphs. We denote by A the adjacency matrix of G. Each node can possibly have a feature vector of size f , which stacks up to an n × f feature matrix X. The degree matrix of A is denoted as D = diag (d 1 , d 2 , · · · , d n ) ∈ R n×n , where d i = vj ∈V A ij . We denote the final node embedding matrix as Z, and evaluate it in both the node clustering and the link prediction tasks. The node clustering task requires the model to partition the nodes into c disjoint groups G 1 , G 2 , · · · , G c , where similar nodes should be in the same group. The target of the link prediction task is to predict whether an edge exists between given node pairs.

Graph Convolutional Network.
Based on the assumption that locally connected nodes are likely to enjoy high similarity (McPherson et al., 2001), each node in most GNN models iteratively smooths the representations of its neighbors for better node embedding. Below is the formula of the l-th GCN layer (Kipf & Welling, 2016a): (1) where A = A+I n is the adjacency matrix of the undirected graph G with self loop added, X (l) is the node embedding matrix at layer l, X (0) is the original feature matrix, Θ (l) are the trainable weights, and δ is the activation function.
A is the smoothing matrix that helps each node to smooth representations of neighboring nodes.
As shown in Eq. 1, each GCN layer contains two operations: feature aggregation (smoothing) and feature transformation. Figure 1 shows the framework of a two-layer GCN. The l-th layer in GCN firstly executes feature smoothing on the node embedding X (l−1) . Then, the smoothed feature X (l−1) is transformed with trainable weights Θ (l) and activation function δ to generate new node embedding X (l) . Note that GCN will degrade to MLP if feature smoothing is removed from each layer.

GNN-based Graph Representation Learning.
GAE (Kipf & Welling, 2016b), the first and the most representative GNN-based graph embedding method, adopts an encoder to generate node embedding matrix Z with inputŝ A and X. A simple inner product decoder is then used to reconstruct the adjacency matrix. The final training loss of GAE is the binary entropy loss between A and A, the reconstructed adjacency matrix and the original adjacency matrix with self loop added. Specifically, the loss function can be defined as follow Where A = sigmoid(Z · Z T ) is the reconstructed adjacency matrix. Motivated by GAE, lots of GNN-based graph representation learning methods have been proposed recently. MGAE (Wang et al., 2017) presents a denoising marginalized autoencoder that reconstructs the node feature matrix X. ARGA (Pan et al., 2018) adopts the adversarial learning strategy, and its generated node embeddings are Figure 1. The framework of a two-layer GCN models.
forced to match a prior distribution. DAEGC (Wang et al., 2019) exploits side information to generate node embeddings in a self-supervised way. AGC (Zhang et al., 2019) proposes an improved filter matrix to better filter out the high-frequency noise. AGE (Cui et al., 2020) further improves AGC by using the similarity of embedding rather than the adjacency matrix to consider original node feature information. Compared with GNN-based graph representation learning methods that rely on the trainable parameters to learn node embeddings, our NAFS is training-free and thus enjoys higher efficiency and scalability.

Observation and Insight
In this section, we make a quantitative analysis on the oversmoothing issue at the node level and then provide some insights when designing NAFS on graphs.

Feature Smoothing in Decoupled GNNs
Recently, many works (Wu et al., 2019;Zhu & Koniusz, 2021;Chen et al., 2020;Zhang et al., 2021c) propose to decouple the feature smoothing and feature transformation in each GCN layer for scalable node classification. Concretely, they execute the feature smoothing operation in advance, and the smoothed features are then fed into a simple MLP to generate the final predicted node labels. The predictive node classification accuracy of these methods is comparable with or even higher than the one of coupled GNNs.
These decoupled GNNs contain two parts: feature smoothing and MLP training. Feature smoothing aims to combine the graph structural information and node features into better features for the subsequent MLP; while MLP training only takes in the smoothed feature and is specially trained for a given task. As stated by previous decoupled GNNs (Wu et al., 2019;Zhu & Koniusz, 2021;Zhang et al., 2021b), the true success of GNNs lies in feature smoothing rather than feature transformation. Correspondingly, we propose to remove feature transformation and preserve the key feature smoothing part alone for simple and scalable node representation.
There is another branch of GNNs that also decouple the feature smoothing and feature transformation. The most representative method of this category is APPNP (Klicpera et al., 2018). It first feeds the raw node features into an MLP to generate intermediate node embeddings. Then the personalized PageRank-based propagation operations are performed on the node embeddings to produce final prediction results. However, compared with scalable decoupled GNNs mentioned in the previous paragraph, this branch of GNNs still have to recursively execute propagation operations in each training epoch, which makes it impossible to perform on large-scale graphs. In the remaining part of this paper, the terminology "decoupled GNNs" refers particularly to the scalable decoupled GNNs mentioned in the previous two paragraphs.

Measuring Smoothing Level
To capture deep graph structural information, a straightforward way is to simply stack multiple GNN layers. However, a large number of feature smoothing operations in a GNN model would lead to indistinguishable node embeddings, i.e., the over-smoothing issue (Li et al., 2018). Concretely, if we executeÂX for infinite times, the node embeddings within the same connected component would reach a stationary state. When adoptingÂ = D r−1Ã D −r ,Â ∞ followŝ which shows that the influence from node v i to v j is only determined by their degrees. Under the extreme condition that r = 0, all the nodes within one connected component have exactly the same representation, making it impossible to apply the node embeddings to subsequent tasks.
Here we introduce a new metric, "Over-smoothing Distance", to measure each node's smoothing level. A smaller value indicates that the node is closer to the stationary state, i.e., closer to over-smoothing.
Definition 1 (Over-smoothing Distance). The Oversmoothing Distance D i (k) parameterized by node i and smoothing step k is defined as where [Â k X] i denotes the i th row ofÂ k X, representing the representations of node v i after smoothing k times; [Â ∞ X] i denotes the i th row ofÂ ∞ X, representing the stationary state of node v i ; Dis(·) is a distance function or a function positively relative with the difference, which can be implemented using Euclidean distance, the inverse of cosine similarity, etc.

Diverse Smoothing Speed across Nodes.
To examine the factors that affect D i (k), we divide nodes in the PubMed dataset into three different groups according to their degrees. In Figure 2, we show the trends of nodes in the first group and the second group where the groupaveraged D i (k) changes with the number of smoothing step increases. The trend of D i (k) averaged over all the nodes is also provided, and the Euclidean distance is chosen as the distance function. Figure 2 shows that the D i (k) of nodes with degrees larger than 60 drops more quickly than nodes with degrees smaller than 3, which implies that nodes with high degrees approach the stationary state more rapidly than the nodes with low degrees.
The quantitative analysis empirically illustrates that the degree of each node plays an essential role in one's optimal smoothing step. Intuitively, nodes with high degrees should have relatively small smoothing steps than nodes with low degrees. In addition, in Appendix A, we have made a detailed theoretical analysis about the graph sparsity, another factor that influences the smoothing speed.

Design Insights of NAFS
Though using feature smoothing operations inside the previous decoupled GNNs is scalable for large graph representation learning, it will lead to sub-optimal node representation. According to the observation in Sec. 3.3, it is sub-optimal to execute feature smoothing for all the nodes indiscriminately as previous decoupled GNNs do since nodes with different structural properties have diverse smoothing speeds. There-fore, node-adaptive feature smoothing (i.e., different nodes have different smoothing levels or mechanisms with equivalent effect) must be adopted to satisfy each node's diverse needs of smoothing level. In our proposed method NAFS, we utilize the metric D i (k) defined in Def. 1 to help accomplish the aim of node-adaptive feature smoothing.

Proposed method
In this section, we present NAFS, a training-free method for scalable graph representation learning. We first compute the smoothed features with the feature smoothing operation. Then the feature ensemble operation is used to combine the smoothed features generated by different smoothing strategies. The pseudo code of NAFS is provided in Appendix D.

Node-Adaptive Feature Smoothing
Figure 3 provides an overview of NAFS. When repeatedly executing X (l) =ÂX (l−1) , the smoothed node embedding matrix X (l−1) contains deeper graph structural information with l increases. The multi-scale node embedding matrices are then combined into a single matrixX such that both local and global neighborhood information are reserved.
The analysis in Sec. 3.3 illustrates that the speed each node achieves its stationary state is highly diverse, which suggests that nodes should be treated individually. To this end, we define "Smoothing Weight" based on D i (k) introduced in Def. 1 for each node so that the smoothing operation can be performed in a node-adaptive manner.
where K is the maximal smoothing step.
To calculate D i (k) more efficiently, an alternative is to replace [Â ∞ X] i in Eq. 4 with X i and implement Dis(·) as the cosine similarity. Larger D i (k) in this case means that node v i is farther from the stationary state and [Â k X] i intuitively contains more relevant node information. Therefore, for node v i , the smoothed feature with larger D i (k) (i.e., larger w i (k)) should contribute more to the final node embedding. The smoothing weight can be formulated in the following matrix form: Weight Matrix W(k) parameterized by smoothing step k is defined as the diagonal matrix derived from η(k) ∈ R n : To achieve the same effect, the feature ensemble operation has multiple knowledge extractors: we vary the value of r in the normalized adjacency matrixÂ r = D r−1 A D −r to easily acquire different knowledge extractors. These knowledge extractors are adopted inside the feature smoothing operation to generate different smoothed features. Concretely, the value of r controls the normalized weight of each edge. So different values of r generate different weight values for all the edges in the graph, which would increase the diversity of our smoothed features evidently. The ablation study in Sec. 6.6 shows that the increased diversity contributes a lot to the high performance of our generated node embeddings when applied to downstream tasks.

Theoretical Analysis
The essential kernel of NAFS is the Smoothing Weight, which determines the output results. We now analyze the factors affecting the value of Smoothing Weight. To simplify our analysis, we suppose r = 0 in the normalized adjacency matrix and apply Euclidean distance as the distance function in Definition 1. Thus we have where || · || 2 symbols two-norm. Theorem 1. For any node i in graph G, there always exists whered i = d i + 1, ||X j || 2 denotes the two-norm of the initial feature of node j, and 0 < λ 2 < 1 denotes the second largest eigenvalue of the normalized adjacency matrixÂ.
Eq. 8 shows that the nodes with smaller degrees may have larger D i (k). Combined with definition 2, we infer that larger D i (k) makes the max k D i (k) more dominant after the softmax operation, causing that weighted average results depend more on itself and its near neighbors. Inversely, for the nodes with smaller degrees, its result of weighted average depends more equally on all itself, its near neighbors, and its distant neighbors.
Besides, D i (k) in a sparser graph (λ 2 is positively relative with the sparsity of a graph) decays slower as k increases. Thus, the weighted average result depends more equally on itself, its near neighbors and its distant neighbors. While for the nodes in a denser graph, the weighted average result depends more on its near neighbors and itself.
The detailed proof of Eq. 8 is in Appendix A.1. Besides, we theoretically analyze how NAFS prevent over-smoothing in Appendix A.2. Lastly, we also theoretically show how NAFS leverages the multiple features over different smooth steps in a node-adaptive manner in Appendix A.3.

Advantages over Traditional Approaches
NAFS generates node embeddings in a training-free manner, making it highly efficient and scalable. Moreover, the node-adaptive smoothing strategy enables it to capture deep structural information. In this subsection, we analyze the advantages of our NAFS over GAE and its variants.
Deep Structural Information. By assigning each node with personalized smoothing weights, NAFS can gather deep structural information without encountering the oversmoothing issue and keep the time and memory cost low. While for GAE and its variants, they either 1) have a coupled structure that cannot go deep due to low efficiency and high memory cost (e.g., GAE) or 2) are unable to adaptively capture structural information and encounter the oversmoothing issue when going deeper (e.g., AGE).
Efficiency. Compared with GAE and its variants, our proposed NAFS does not have any trainable parameters, giving it a significant advantage in efficiency. When generating node embeddings, NAFS only needs to execute the feature smoothing and feature ensemble operations, which has a time complexity of O(Kmf + Kn), where K denotes the maximal smoothing step. On the other hand, the asymptotic training time complexity of GAE and its variants is O(nf 2 + n 2 f ), which is one magnitude larger than NAFS. In Sec. 6.4, we further present efficiency comparison in detail on real-world datasets between NAFS and GAE.
Memory Cost. Our proposed NAFS enjoys not only high efficiency but also low memory cost. NAFS only requires to store the sparse adjacency matrix and the smoothed features, and thus the memory cost is O(m + nf ), which grows linearly with graph size n in typical real-world graphs. For GAE and its variants, they need additional space to store the training parameters and especially the reconstructed dense adjacency matrix, making their memory cost O(n 2 + m + nf ), which is also one magnitude larger than NAFS. Scalability comparison conducted on synthetic datasets between NAFS and GAE can be found in Appendix C.2.

Experimental Results
In this section, we conduct extensive experiments to evaluate the proposed NAFS. We first introduce the considered baseline methods, used datasets, and experimental settings. Then, we demonstrate the advantages of NAFS from the following three perspectives: (1) end-to-end comparison with the state-of-the-art methods, (2) scalability and efficiency, and (3) effectiveness along with ablation studies.

Datasets and Baseline Methods
Several widely-used network datasets (i.e., Cora, Citeseer, PubMed, Wiki and ogbn-products) are used in our experiments. We include the properties of these datasets in Appendix B.1.
For different tasks, the compared baselines are as follows: •   For NAFS, we investigate the following three variants: NAFS-mean, NAFS-max, and NAFS-concat. They all have the complete NAFS framework and only differ in the ensemble strategy adopted in feature ensemble operation.

Experimental Settings
Node Clustering. For the node clustering task, we apply K-Means (Hartigan & Wong, 1979) to node embeddings to get the clustering results. Three widely-used metrics are used for evaluation: Accuracy (ACC), Normalized Mutual Information (NMI), and Adjusted Rand Index (ARI).
Link Prediction. For the link prediction task, 5% and 10% edges are randomly reserved for the validation set and the test set. Once the node embeddings have been generated, we follow the same procedure as GAE to reconstruct the adjacency matrix. Two metrics -Area Under Curve (AUC) and Average Precision (AP) are used in the evaluation of the link prediction task.
We run the compared baselines with 200 epochs and repeat the experiment 10 times on all the datasets, and report the mean value of each evaluation metric. The detailed setting of the hyperparameters and experimental environment are introduced in Appendix B.2 and B.3.

End-to-end Comparison
Link Prediction. Table 1 shows the performance of different methods on the link prediction task. On the three datasets, the proposed NAFS consistently achieves the best or the second-best performance compared with all the baseline methods. Remarkably, NAFS-concat achieves stateof-the-art performance on PubMed, and outperforms the current state-of-the-art method -ARGA by 0.8% and 0.1% on AUC and AP, respectively.
Node Clustering. The node clustering results of each method are shown in Table 2. Among three ensemble strategies, NAFS-concat has the overall best performance across the three datasets, which also consistently outperforms the strongest baseline -AGE. For example, the best of NAFS variants exceeds AGE by 2.6%, 1.8%, 0.7%, and 3.5% on   Cora, Citeseer, PubMed, and Wiki, respectively. Although NAFS-concat has the overall best performance, it falls behind NAFS-mean or NAFS-max on some metrics of some datasets; and it needs more CPU memory to store the node embeddings and longer time for inference.
It is quite interesting to find that the three NAFS variants outperform training-based GAE on all the four datasets, and they outperforms the current SOTA method, AGE, in most cases. The competitive performance of NAFS further illustrates that it is possible to achieve decent performance on graphs with only feature smoothing without any training parameters.

Efficiency Analysis
To validate the efficiency of NAFS, we compare it with GAE and AGE, measuring their overall running time for generating node embeddings. The comparison is conducted on the three citation networks. For fairness, all methods only use CPUs for computation. Figure 5 showcases the results along with the speedup ratio of NAFS against AGE. Figure 5 shows that NAFS is significantly faster than the considered methods on the three datasets. For example, NAFS is almost two magnitudes faster than AGE on the relatively large dataset -PubMed. The high efficiency of NAFS is attributed to no trainable parameters, while GAE and AGE are both training-based methods. Although AGE executes the feature smoothing step in advance, it adopts a time-consuming ranking policy to select positive and negative examples. This alteration in AGE brings both superior performance and longer running time compared with GAE.

Scalability Analysis
A challenging task in graph representation learning is how to effectively get decent node embedding on very large graphs. To verify the scalability advantage of NAFS, we add the evaluation on the large real-world datasets. The experiment settings are altered compared to Sec. 6.4: both GAE and AGE are trained on an NVIDIA TITAN RTX, which has 24 GB of memory since only using CPU to train is unacceptable on large graphs in real-world applications.
Concretely, we sample subgraphs of different scales from the full ogbn-products graphs via uniformly random node selection and then report the sampled graph size (i.e., number of nodes) and the corresponding runtime (seconds) in Table 3. The experimental results show that NAFS can support the larger graphs (i.e., larger than 30,000 nodes) than the compared baselines. Besides, it is significantly faster than the compared baselines, especially for large graph datasets. We also include more details in Appendix C.2 to demonstrate the excellent scalability of NAFS (See Figure 7).

Ablation Study
To thoroughly investigate the proposed NAFS, ablation studies on the node clustering task are designed to analyze the effectiveness of feature smoothing and feature ensemble in NAFS. The experiments are conducted on the PubMed dataset, and Normalized Mutual Information (NMI) is used to measure the performance.
Different Weighting Strategies. One important operation in NAFS is to use the node-adaptive weight to average the smoothed features of different smoothing steps. Here we change the "Adaptive Weight" in our method to "Single Hop" (only the smoothed feature at the last smoothing step is reserved) and "Naive Average" (smoothed features at every smoothing step has equal weight) and evaluate their performance. Figure 4(a) shows the performance of these three different weighting strategies.
From Figure 4(a), the weighting strategy "Single Hop" performs the worst among the three, since it only makes use of the information at the last smoothing step, which could lead to the over-smoothing issue when the maximal smoothing step becomes large. On the other hand, the weighting strategy "Adaptive Weight" shows a better performance than "Naive Average". Besides, when the number of maximal smoothing steps becomes large, the performance of "Naive Average" begins to drop, while "Adaptive Weight" does not. "Naive Average" assigns the same weight across all the different steps of smoothed features, which also leads to over-smoothing when it tries to exploit extremely deep structural information. Instead, by assigning adaptive weights, the "Adaptive Weight" strategy in NAFS could exploit such deep information and avoid the over-smoothing issue, thus improving the quality of the generated node embedding.
Different Ensemble Strategies. We use different ensemble strategies to obtain the final node embeddings in our proposed method -NAFS, which include "Mean", "Max", and "Concat". To evaluate the impacts of these different ensemble strategies, we change the maximal smoothing step, K, from 15 to 40, and evaluate corresponding node clustering performance. For reference, the performance of NAFS without feature ensemble, "r=0.3 only" is also reported. The experimental results in Figure 4(b) illustrate that at most times, "r=0.3 only" performs worse than others, which shows the effectiveness of the three different ensemble strategies. However, if we limit the comparison within the ensemble strategies, the performance superiority is unclear. In Table 1 and 2, we can also have that different ensemble strategies perform diversely across different datasets and tasks. It indicates that these three ensemble strategies all have necessities in their own way.

Interpretability
To better understand why NAFS is effective, we visualize the node embedding generated by NAFS using T-SNE (Van der Maaten & Hinton, 2008) on the Cora dataset. Moreover, we also visualize the node embedding ofÂ k X for reference. All the visualization results are shown in Figure 6. The first row is the results of our proposed NAFS, and the second row is the results ofÂ k X.
Figure 6(c) and 6(k) illustrate that at hop 10, the node embedding produced by NAFS andÂ k X are both distinguishable. But as the value of maximal smoothing step becomes larger, the node embedding ofÂ k X falls into total disorder, like the situation showed by Figure. 6(n), 6(o) and 6(p). At the same time, NAFS is able to maintain the distinguishable results even when the value of the maximal number of smoothing step increases to 200.

A.1. What influences the Smoothing Weight?
The essential kernel of NAFS is the Smoothing Weight, which determines the output results. We now analyze the factors affecting the value of Smoothing Weight. To simplify our analysis, we suppose r = 0 in the normalized adjacency matrix and apply Euclidean distance as the distance function in Definition 1. Thus we have where || · || 2 symbols two-norm.
Theorem 2. For any node i in graph G, there always exists whered i = d i +1,d j = d j +1, ||X j || 2 denotes the two-norm of the initial feature of node j, and 0 < λ 2 < 1 denotes the second largest eigenvalue of the normalized adjacency matrixÂ.
To prove Theorem 2, we introduce the following lemma. Lemma 3.
where ei denotes a one-hot row vector with its i th components as 1 and other components as 0, λ2 is the second largest eigenvalue ofÂ anddi denotes the degree of node i plus 1 (to include itself).
The proof of Lemma 3 can be found in (Chung & Graham, 1997). Next we will prove Theorem 2.
Proof. According to equation 7, we can have that where X jp denotes the p th feature of node j.
Based on Lemma 3 we then consider the smoothing distance for weighted averaged embedding features to the station state.
Therefore there holds Theorem 2 We denote the constant n j=1 (d j ||X j || 2 2 ) as cdx because it is independent with j, then Theorem 2 can be written as We then analyze the factors affecting the Smoothing Weight on a specific node v i . From Eq. 15 we know that the nodes with smaller degrees may have larger D i (k). Combined with definition 2, we infer that larger D i (k) makes the max k D i (k) more dominant after the softmax operation, causing that weighted average results depend more on itself and its near neighbors. Inversely, for the nodes with smaller degrees, its result of weighted average depends more equally on all itself, its near neighbors, and its distant neighbors.
At the same time, the Smoothing Weight of the node in a sparser graph (λ 2 is positively relative with the sparsity of a graph) decays slower as k increases. Thus, the weighted average result depends more equally on itself, its near neighbors and its distant neighbors. While for the nodes in a denser graph, the weighted average result depends more on its near neighbors and itself.

A.2. How NAFS prevent over-smoothing?
Based on Theorem 2 we then consider the smoothing distance for weighted averaged embedding features to the station state: When ω i (k) = 1/(K + 1), which means the average option is non-weighted, we have causing over-smoothing. When In real-world graphs, nodes have different initial features, thus there is little chance that the combination of a node's neighboring features both lies in the opposite direction and is of the same norm of the node's initial feature. Suppose that there exits a constant > 0 satisfying Thus we see that even K goes to infinity, we are able to prevent the node representations from reaching the stationary state (the distance bound depends on the node degree d i and the initial feature X i ). Note that in practice the representation at k-th smoothing step X (k) achieves the stationary state much earlier than the infinity-th smoothing step we use in the theoretical analysis, so we use the softmax normalization in Equation 5 to produce a slightly larger bias towards features with longer distances to the stationary state.
The above analysis theoretically proved our weighting scheme is able to prevent over-smoothing as the smoothing step goes to infinity. We now show that another advantage of our method is that it can fully leverage the multiple features over different smooth steps in a node-adaptive manner, which is different from the traditional routine of the smoothblind and fixed scheme. Suppose that the feature of k step is not over-smoothed yet for a specific node i, i.e., the distance of D i (k) = i > 0, we have: We see that as long as the feature is not over-smoothed, our method will assign a non-zero weight to the feature. Further, we see that the bound of weight can be affected by the oversmoothing distance of i and degree d i of a specific node, implying an adaptive weighting strategy.

B.1. Datasets Description
Cora, Citeseer, PubMed, and Wiki are four popular network datasets. The first three (Yang et al., 2016) are citation networks where nodes stand for research papers, and an edge exists between a node pair if one cites the other. Wiki (Yang et al., 2015) is a webpage network where nodes stand for webpages, and an edge exists between a node pair if one links the other. The ogbn-arxiv (Hu et al., 2020) dataset is also a citation network that contains more than 160M nodes. Besides, the ogbn-products dataset is an undirected and unweighted graph, representing an Amazon product co-purchasing network. Table 4 presents an overview of these six datasets.

B.2. Hyperparameters setting
When generating node embeddings, we use the values of r in [0, 0.1, 0.2, 0.3, 0.4, 0.5] to get six different normalized adjacency matrixÂ. The only exception is the link prediction task on the PubMed dataset, where we use values of r in [0.3, 0.4, 0.5] instead. The optimal value of the maximal smoothing steps ranges from 1 to 70. Hyperparameters for all the baseline methods are tuned with OpenBox (Li et al., 2021) or following the settings in their original paper.

B.3. Experimental Environment
The experiments are conducted on a machine with Intel(R) Xeon(R) Gold 5120 CPU @ 2.20GHz, and a single NVIDIA TITAN RTX GPU with 24GB GPU memory. The operating system of the machine is Ubuntu 16.04. For software versions, we use Python 3.6, Pytorch 1.7.1, and CUDA 10.1.

C.1. Performance on the Node Classification Task
Node classification performance. In this part, we assess the quality of the node embeddings generated by NAFS by the evaluation on the node classification task. We follow the linear evaluation protocol, which applies a linear classifier (i.e., Logistic Regression) to the node embeddings to generate final prediction results. We choose GCN (Kipf & Welling, 2016a), JK-Net (Xu et al., 2018), C&S (Huang et al., 2020), SGC (Wu et al., 2019), GAT (Veličković et al., 2017), PPRGo (Bojchevski et al., 2020), APPNP (Klicpera et al., 2018), andDAGNN (Liu et al., 2020) as comparison baselines on the node classification task. Note that C&S in our evaluation adopts a two-layer MLP as the base model. The evaluation results on three popular datasets, Cora, Citeseer, and PubMed (Yang et al., 2016), are provided in Table 5.
The table shows that our method achieves comparable predictive accuracy as one of the state-of-the-art methods DAGNN. However, since the node embeddings are generated before training, our method avoids performing recursive feature smoothing at each training epoch and storing the entire adjacency matrix on GPU. In this way, our method is more scalable and efficient to apply on large graphs than most state-of-the-art methods like DAGNN.
Performance-efficiency comparison. In this part, we evaluate the node classification accuracy of NAFS-mean and several baseline methods on the ogbn-arxiv dataset (Hu et al., 2020). The predictive accuracy and runtime results are both shown in Figure 8 whose X-axis is in log scale. The Figure shows that NAFS-mean only falls behind GAT and outperforms many competitive baselines in terms of test accuracy. Although the test accuracy of GAT is slightly higher than NAFS-mean, NAFS-mean achieves over 6x speedup than GAT.

C.2. Scalability Comparison on Synthetic Graphs
To test the scalability of our proposed NAFS, we also use the Erdős-Rényi graph generator in the Python package NetworkX (Hagberg et al., 2008) to generate artificial graphs of different sizes. The node sizes of the generated artificial graphs vary from 5,000 to 100,000, and the probability of an edge exists between two nodes is set to 0.0001. The experiment settings are altered compared to Sec. 6.4: both GAE and AGE are trained on an NVIDIA TITAN RTX, which has 24 GB of memory since only using CPU to train is unacceptable on large graphs in real-world applications.
The overall experiment results are shown in Figure 7. Figure 7(a) shows a more detailed comparison on relatively small artificial graphs whose node sizes range from 5,000 to 20,000. Besides, Figure 7(b) shows that GAE encounters the out-of-memory problem on the artificial graph composed of 30,000 nodes, while AGE encounters the out-of-memory problem at graph size 80,000. Our proposed NAFS is successfully carried out on all the artificial graphs. On the artificial graph consisting of 100,000 nodes, NAFS accomplishes the node embedding generation in 66.1 seconds, which is less than the running time of AGE on the artificial graph consisting of 10,000 nodes, 80.9 seconds. The graph size limit of our proposed NAFS is bounded by the CPU memory size. As long as one can successfully execute the multiplication of the sparse adjacency matrix and the feature matrix, NAFS can then be implemented.

C.3. Performance on sparse graphs
To validate the performance benefits of NAFS on sparse graphs, we conduct experiments on the PubMed dataset under two handcrafted sparsity settings: edge sparsity and feature sparsity.
Two sparsity settings. Under the edge sparsity setting, we randomly remove some edges in the original PubMed dataset to strengthen the edge sparsity issue. The edges removed from the original graph are kept the same across all the compared methods under the same edge removing rate. Under the feature sparsity setting, we randomly choose some nodes in the original PubMed dataset and set their feature vectors to all-zero vectors. The selected nodes are kept the same across all the compared methods under the same feature removing rate.  Experiment settings. We report the performance of GAE and NAFS-mean on the link prediction task under different sparsity settings. We adopt two metrics -Area Under Curve (AUC) and Average Precision (AP) to evaluate the performance of each method. The evaluations are conducted with the edge removing rate and the feature removing rate set to 0.2, 0.4, 0.6, respectively. We repeat each method 10 times and report the mean performance and the corresponding standard deviations in Table 6 and 7.
Experiment results. The experimental results from Table 6 and 7 show that 1) under both edge and feature sparsity settings, NAFS consistently outperforms GAE and GCN on the link prediction task and the node classification task, respectively. 2) the performance gains are larger in sparser graphs with larger feature/edge removing rates. Concretely, the AUC of NAFS outperforms GAE by a margin of 1.0% if the edge is not removed, and the performance gain has increased to 3.7% when the edge removing rate is 0.6.  The experiment results illustrate that NAFS can effectively exploit distant neighborhood information without the oversmoothing problem and thus get better performance on sparse graphs compared with baseline methods.

C.4. Comparison with NDLS
NAFS differs from the related work NDLS in: 1) NDLS searches for the optimal smoothing step for each node, while NAFS presents "smoothing weight" to aggregate the outputs after different smoothing steps. 2) Compared with NDLS, NAFS further proposes to combine the smoothed outputs of different smoothing operators. 3) NDLS only focuses on the semi-supervised node classification task, while NAFS further considers the unsupervised task. The comparison in Table 8 shows NAFS's performance superiority over NDLS on unsupervised tasks.

D. More details of NAFS D.1. Pseudo code of NAFS
Alg. 1 shows the whole pipeline of our proposed NAFS. We first initialize X (0) as the original feature matrix X. Given the normalization parameter r t , we obtained the corresponding normalized adjacency matrixÂ rt = D rt−1 A D −rt , which acts as knowledge extractors (line 4). After that, for each node, we use E.q. 4 to calculate the Over-smoothing Distance with all the k ranging from 0 to K (line 6, 7). Then we calculate its Aggregation Weights through Eq. 5 (line 8, 9). After obtaining Aggregation Weights with all k and i, we construct the Aggregation Weight matrix for each k with Eq. 6 (line 10, 11). Next, we compute the NAFS out-putX (t) withX (t) = K k=0 W(k)Â k rt X (line 12). Finally, we compute the final embedding result of all t through Z ← ⊕ i∈{1,2,...,T }X (i) (line 14).

D.2. Motivation of Feature Ensemble
Adopting different smoothing operators in the feature smoothing operation (the normalized adjacency matrixÂ in Eq. 1) is equivalent to smoothing features in different manners. Other than the normalized adjacency matrixÂ, there are many alternatives that have been proposed recently. For example, GraphSAGE (Hamilton et al., 2017) designs three smoothing operators (i.e., Mean, LSTM and Pooling) to flexibly capture the information of neighboring nodes. SIGN (Frasca et al., 2020) enriches the smoothing operators with Personalized-PageRank-based (Klicpera et al., 2019) and triangle-based (Monti et al., 2018) adjacency matrices. However, these methods are not designed for graph representation learning, and different smoothing operators may result in diverse smoothed features. For better node representation, the feature ensemble is used to combine the smoothed feature under different smoothing operators.

D.3. Relations with other Scalable GNN Architectures
The scalable GNNs can be roughly classified into two categories: (a) "first propagate then predict"; (b) sampling + ordinary GNN. The first category includes famous GNN models like SGC (Wu et al., 2019) and SIGN (Frasca et al., 2020). They disentangle the coupled propagation and transformation operations in traditional GCN layers (Kipf & Welling, 2016a), and execute all the propagation operations as preprocessing. The most representative GNN model of the second category is GraphSAGE (Hamilton et al., 2017), and many other works have been proposed following it, like FastGCN (Chen et al., 2018) and GraphSAINT (Zeng et al., 2019). The main contribution of these works is effective sampling strategies that can preserve the most valuable information from the original graph. The sampling strategies always act as a plug-and-play module, and can be combined with ordinary GNNs, which allows the latter to perform on large-scale graphs. GNNs belong to the "first propagate then predict" category usually enjoy higher efficiency since they avoid recursively performing propagation in each training epoch. Our proposed NAFS belongs to the "first propagate then predict" category.

D.4. Advantage in Distributed Settings
NAFS can also be adapted to the distributed environment. The feature smoothing process in NAFS is sparse matrix dense matrix multiplications, which have mature implementations in distributed environments. Besides, this process only needs to be pre-computed at once. In contrast, during each training iteration of GAE and its variants, each node must repeatedly pull the intermediate representations of other nodes, leading to high communication costs.

D.5. Association with Mixing Time in Markov Chain
Setting r = 0 in the adjacency matrix makes the feature smoothing process a Markov chain. To study the mixing time at the node level, we denote t mix ( , i) as the minimum step for node i such that the distance between its representation and its stationary state is at most . It can be calculated as follows with the help of Lemma 3 in Appendix A.1: