How Graph Neural Networks Learn: Lessons from Training Dynamics
A long-standing goal in deep learning has been to characterize the learning
behavior of black-box models in a more interpretable manner. For graph neural
networks (GNNs), considerable advances have been made in formalizing what
functions they can represent, but whether GNNs will learn desired functions
during the optimization process remains less clear. To fill this gap, we study
their training dynamics in function space. In particular, we find that the
gradient descent optimization of GNNs implicitly leverages the graph structure
to update the learned function, as can be quantified by a phenomenon which we
call \emph{kernel-graph alignment}. We provide theoretical explanations for the
emergence of this phenomenon in the overparameterized regime and empirically
validate it on real-world GNNs. This finding offers new interpretable insights
into when and why the learned GNN functions generalize, highlighting their
limitations in heterophilic graphs. Practically, we propose a parameter-free
algorithm that directly uses a sparse matrix (i.e. graph adjacency) to update
the learned function. We demonstrate that this embarrassingly simple approach
can be as effective as GNNs while being orders-of-magnitude faster.
Chenxiao Yang, Qitian Wu, David Wipf, Ruoyu Sun, Junchi Yan
# How Graph Neural Networks Learn: Lessons from Training Dynamics in Function Space

###### Abstract

A long-standing goal in deep learning has been to characterize the learning behavior of black-box models in a more interpretable manner. For graph neural networks (GNNs), considerable advances have been made in formalizing what functions they can represent, however it remains less clear whether and how GNNs learn desired functions during the optimization process. To fill this critical gap, we study the learning dynamics of GNNs in function space via the analytic framework of overparameterization. In particular, we find that the seemingly complicated training process of GNNs can be re-cast into a more familiar label propagation framework, due to the graph inductive bias implicit in this process. From this vantage point, we provide explanations for why the learned GNN functions successfully generalize and for their pathological behavior on heterophilic graphs, which are consistent with observations. Practically, sparsifying and implementing the learning dynamics lead to a minimalist semi-supervised learning algorithm with the efficiency of classic algorithms and the effectiveness of modern GNNs.

## 1 Introduction

_Graph Neural Networks (GNNs)_(Gori et al., 2005; Scarselli et al., 2008; Bruna et al., 2014; Kipf and Welling, 2017) represent network architectures for learning on entities with explicit relations. In addition to their empirical success, the recent pursuit of theoretical understanding has also led researchers to dissect GNN models, especially regarding their representation (Maron et al., 2019; Xu et al., 2019; Oono and Suzuki, 2019; Chen et al., 2019) and generalization (Scarselli et al., 2018; Verma and Zhang, 2019; Garg et al., 2020) powers. Another research area in parallel with representation and generalization is optimization, which focuses on understanding how the training algorithm identifies desirable GNN functions. Despite its importance, the optimization process of training GNNs is far less understood compared with other topics. While previous work has studied the learning dynamics of linearized GNNs (Xu et al., 2021) in their _weight space_, with emphasis on the convergence speed of empirical risk, an analysis of their dynamics in _function space_ that directly examines how the GNN function evolves has been lacking. In particular, we focus on the following open questions:

_1 How do GNN models evolve in function space during optimization, and how exactly do graphs regularize or bias this process?_

_2 In light of the above, are there interpretable ways to explain why the learned GNN function successfully generalizes and when they might fail?_

Studying these questions can lead to deeper understanding of GNNs and provide guidelines for designing more principled algorithms (as we will exemplify in Sec. 3). As an initial probe in this direction, we study GNN learning dynamics in function space, i.e., how the learned function \(f_{t}(\mathbf{x})\) (or model prediction) evolves on an arbitrary sample \(\mathbf{x}\) as model weights \(\mathbf{W}_{t}\) evolves over time according to gradient descent. Our contributions are summarized as follows:
**Theoretical Contributions.** We provide answers to the above two questions, via the widely adopted analytical framework (Jacot et al., 2018; Arora et al., 2019; Du et al., 2019) of studying neural networks in overparameterized regimes, and map our insights to practice via empirical verification.

\(\circ\)_GNN Dynamics Resemble Label Propagation:_ The evolution of GNNs in function space follows a cross-instance propagation scheme with pairwise weights controlled by a kernel matrix \(\mathbf{\Theta}_{t}\) called the _Neural Tangent Kernel (NTK)_, which we demonstrate tends to align well with, or even become equivalent to, various specialized graph adjacency matrices \(\mathbf{A}\) in overparameterized regimes. In this way, each gradient descent step for GNNs resembles classic _Label Propagation (LP)_ algorithms (Szummer and Jaakkola, 2001; Zhu and Ghahramani, 2002; Zhu et al., 2003; Zhou et al., 2003; Chapelle et al., 2009) that propagate true labels from the training set to unseen samples along the graph, as illustrated in Fig. 1(a). Such a characterization of the evolution of learned GNN functions provides answers to question 1 and reveals the role of graphs in biasing this procedure.

\(\circ\)_Rethinking Generalization and Heterophily:_ The alignment between NTK \(\mathbf{\Theta}_{t}\) and graph adjacency \(\mathbf{A}\) (dubbed as kernel-graph alignment) explains success and pathology of learned GNN functions. Let \(\mathbf{\Theta}^{*}\) be the optimal kernel denoting if samples share the same (true) label, and the alignment between \(\mathbf{\Theta}^{*}\) and \(\mathbf{A}\) indicates the homophily level of a graph (i.e. whether connected nodes indeed have the same label). As illustrated by the triangle relation in Fig. 1(b), a larger homophily level adding inherently good kernel-graph alignment translates to better kernel-target alignment (Cristianini et al., 2001), which is desired for favorable generalization. This answers question 2 and sheds additional insights on heterophily problems. As support, we establish a strong correlation between generalization and homophily level in overparameterized regimes, and characterize the Bayesian optimality of GNNs that minimizes the population risk. Our empirical results on the time evolution of real-world GNN NTKs on various datasets serves to corroborate our analysis.

**Practical Contributions.** Drawing upon the novel message passing perspective on the learning dynamics of neural networks, we propose to directly replace the NTK matrix with the adjacency matrix \(\mathbf{A}\), which leads to a parameter-free algorithm which we refer to as _Residual Propagation (RP)_. Theoretically, it serves as an extreme-case characterization of the learning dynamics of GNNs where graph structure dominates the optimization process. Empirically, we found this embarrassingly simple algorithm can outperform GNNs by \(3.0\%\) on OGB benchmarks even without using node features (improvement also shown on \(12\) additional datasets in Appendix). And notably RP takes orders of magnitude (\(<0.01\)) less time to excel the performance of GNNs.

## 2 Preliminaries

**Notations and Setup.** Given a set of training inputs \(\mathcal{X}=\{\mathbf{x}_{i}\}_{i=1}^{n_{l}}\in\mathbb{R}^{n_{l}\times d}\) and labels \(\mathcal{Y}=\{y_{i}\}_{i=1}^{n_{l}}\in\mathbb{R}^{n_{l}}\), where \(n_{l}\) is the size of the labeled instances, we aim to learn a predictive function \(f(\mathbf{x})\) (such as a neural network) parameterized by weights \(\mathbf{W}\). For (semi-)supervised learning, we minimize the squared loss \(\mathcal{L}\) using _Gradient Descent (GD)_,

Figure 1: **(a) Learning dynamics of (graph) neural networks from a function space viewpoint where residuals (i.e. difference between ground-truth labels and predictions) propagate from observed to unobserved samples based on a kernel matrix and shape the learned function. (Each image is a data point) (b) Alignment of matrices that control the residual propagation process: NTK \(\mathbf{\Theta}_{t}\), adjacency matrix \(\mathbf{A}\), and the optimal kernel matrix \(\mathbf{\Theta}^{*}\). Kernel-graph alignment is introduced in this work.**
\[\mathcal{L}=\|\mathcal{F}_{t}-\mathcal{Y}\|^{2}/2,\quad\partial\mathbf{W}_{t}/ \partial t=-\eta\nabla\mathbf{W}\mathcal{L}, \tag{1}\]

where \(\mathbf{W}_{t}\) and \(\mathcal{F}_{t}=\{f_{t}(\mathbf{x})\}_{\mathbf{x}\in\mathcal{X}}\in\mathbb{R}^{n_{l}}\) are weights and predictions for the training set at optimization time index \(t\), and \(\eta\) is the learning rate. Temporal discretization of this gradient flow system with time-step \(\Delta t=1\) yields the fixed step-size GD algorithm commonly used in practice, i.e. \(\mathbf{W}_{t+1}=\mathbf{W}_{t}-\eta\nabla\mathbf{W}\mathcal{L}\). Additionally, \(\mathcal{X}^{\prime}\) and \(\mathcal{Y}^{\prime}\) denote testing instances, such that \(\bar{\mathcal{X}}=[\mathcal{X},\mathcal{X}^{\prime}]\in\mathbb{R}^{n\times d}\) (resp. \(\bar{\mathcal{Y}}\)) represents the concatenation of training and testing inputs (resp. labels), where \(n\) is the full dataset size. For convenience, we generally refer to \(f_{t}(\mathbf{x})\) as prediction for a single data point, which is allowed to use additional information such as input features of other nodes by slightly abusing notation, \(\mathcal{F}_{t}\) and \(\mathcal{F}^{\prime}_{t}\) as predictions for the training and testing sets, which are agnostic to model parameterization in order to apply to parameter-free models. Our analysis generalizes to other loss functions and multi-dimensional output (see Appendix. C.2 and C.3).

Following (Xu et al., 2021), we focus on node classification or regression tasks, where instances (i.e. nodes) and their relations (i.e. edges) are described by an undirected graph \(\mathcal{G}=(\mathcal{V},\mathcal{E})\), \(|\mathcal{V}|=n\). The graph defines a symmetric adjacency matrix \(\mathbf{A}\in\mathbb{R}^{n\times n}\) where \(\mathbf{A}_{ij}=\mathbf{A}_{\mathbf{x};\mathbf{x}_{j}}=1\) for a pair of connected nodes \((\mathbf{x}_{i},\mathbf{x}_{j})\) otherwise \(0\). Based on the data split, we denote submatrices of \(\mathbf{A}\) using \(\mathbf{A}_{\mathcal{X}\mathcal{X}}\) and \(\mathbf{A}_{\mathcal{X}^{\prime}\mathcal{X}}\). As a standard technique in graph theory (Chung, 1997), we let \(\mathbf{A}\) be normalized, i.e. \(\mathbf{A}\leftarrow\mathbf{D}^{-\frac{1}{2}}(\mathbf{A}+\mathbf{I})\mathbf{D }^{-\frac{1}{2}}\) where \(\mathbf{D}\) is the node degree matrix. Our insights apply to both transductive and inductive settings (see Appendix. C.1), but will focus on the former case unless stated otherwise.

**Graph Neural Networks (GNNs)** are a class of network architectures for learning representations on graphs, commonly trained based on the aforementioned optimization procedure. For GCN (Kipf & Welling, 2017) (and other GNNs that differ by the definition of \(\mathbf{A}\)), each layer can be written as

\[\mathbf{Z}^{(\ell)}=\sigma(\mathbf{A}\mathbf{Z}^{(\ell-1)}\mathbf{W}^{(\ell) })\in\mathbb{R}^{n\times m}, \tag{2}\]

where \(\mathbf{Z}^{(\ell)}\) are node representations at the \(\ell\)-th layer with \(\mathbf{Z}^{(0)}=\bar{\mathcal{X}}\), \(\sigma\) is ReLU activation, and \(m\) is the model width. We denote the GNN prediction for a single data point as \(f(\mathbf{x};\mathbf{A})\), conditioned on \(\mathbf{A}\) to distinguish it from other models. We adopt this architecture as the default GNN model.

**Label Propagation (LP)** represents a class of algorithms for semi-supervised learning, ground-truth labels \(\mathcal{Y}\) are propagated along edges to predict \(\mathcal{F}^{\prime}\) for the testing (unlabeled) set. From (Zhu & Ghahramani, 2002), the LP update equation can be written as

\[\mathrm{LP}(\mathcal{Y};k,\alpha)=[\mathcal{F}_{k},\mathcal{F}^{\prime}_{k}] =\alpha\mathbf{A}\;[\mathcal{F}_{k-1},\mathcal{F}^{\prime}_{k-1}]+(1-\alpha)[ \mathcal{Y},\mathbf{0}], \tag{3}\]

where \([\cdot,\cdot]\) is concatenation, \(k\) is the iteration number, and \(\alpha\) is a hyperparameter. As initialization \([\mathcal{F}_{0},\mathcal{F}^{\prime}_{0}]=[\mathcal{Y},\mathbf{0}]\), and after convergence \([\mathcal{F}_{\infty},\mathcal{F}^{\prime}_{\infty}]\propto(\mathbf{I}_{n}- \alpha\mathbf{A})^{-1}[\mathcal{Y},\mathbf{0}]\). LP algorithms have found wide applicability due to their superior efficiency, scalability, and simplicity for deployment.

## 3 A Label Propagation View of Gradient Descent Dynamics

Before delving into GNNs specifically, we detour by providing a LP perspective on the evolution of general parameterized models during GD-based optimization. From this vantage point we propose a simple and interpretable propagation algorithm for semi-supervised learning, with an aim towards illustrating core ideas that will later contribute to our understanding of GNN behavior during optimization (Sec. 3.1). This new algorithm has interesting connections with classic LP algorithms, encompassing them as special cases with additional global optimization guarantees. We evaluate on OGB benchmarks (and additional \(12\) datasets in Appendix E) to verify effectiveness, efficiency, and scalability (Sec. 3.2). We will circle back to connections with GNN learning dynamics in Section 4.

### Learning Dynamics in Function Space and Residual Propagation

We start by characterizing the evolution of a general model (with no restriction on model architecture) in its function space \(\mathcal{F}_{t}\), induced by GD-based optimization that continuously updates the weights \(\mathbf{W}_{t}\). On the training set, this process can be described by (Jacot et al., 2018)

\[\partial\mathcal{F}_{t}/\partial t=\eta\;\mathbf{\Theta}_{t}(\mathcal{X}, \mathcal{X})\mathcal{R}_{t},\quad\text{(NTK Matrix): }\;\mathbf{\Theta}_{t}(\mathcal{X},\mathcal{X})\triangleq\nabla\mathbf{w} \mathcal{F}_{t}^{\top}\nabla\mathbf{w}\mathcal{F}_{t}\in\mathbb{R}^{n_{l} \times n_{l}}, \tag{4}\]

where \(\mathcal{R}_{t}=\mathcal{Y}-\mathcal{F}_{t}\in\mathbb{R}^{n_{l}}\) denotes _residuals_(Hastie et al., 2009) (a.k.a. errors), the difference between ground-truth labels and model predictions. The so-called _Neural Tangent Kernel (NTK)_(Jacot
et al., 2018)\(\mathbf{\Theta}_{t}(\mathcal{X},\mathcal{X})\) is produced by the product between Jacobian matrices, which is dependent on the network architecture, and evolves over time due to its association with time-varying weights \(\mathbf{W}_{t}\). Specially, if the kernel is constant (such as for linear models and infinitely-wide neural networks), (4) reduces to the learning dynamics of kernel regression (Shawe-Taylor & Cristianini, 2004).

**Information Propagation Interpretation.** While (4) has found usage for analyzing the convergence of empirical risk (Du et al., 2019; Arora et al., 2019) and the spectral bias of deep learning (Mei et al., 2019; Cao et al., 2019), it is restricted to a limited set of samples (i.e. training set). To see how the model evolves on arbitrary inputs for fully characterizing the learned function, we extend (4) to accommodate unseen samples (which could be chosen arbitrarily). Specifically, let \(\mathcal{R}^{\prime}_{t}=\mathcal{Y}^{\prime}-\mathcal{F}^{\prime}_{t}\) denote residuals for the testing set, and after temporal discretization, (4) can be extended and rewritten neatly using a single variable residual \(\mathcal{R}\) (see derivation in Appendix B.1)

\[\left[\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\right]=-\eta\ \left[ \begin{array}{cc}\mathbf{\Theta}_{t}(\mathcal{X},\mathcal{X})&\mathbf{\Theta}_{t}( \mathcal{X},\mathcal{X}^{\prime})\\ \mathbf{\Theta}_{t}(\mathcal{X}^{\prime},\mathcal{X})&\mathbf{\Theta}_{t}(\mathcal{X }^{\prime},\mathcal{X}^{\prime})\end{array}\right]\left[\mathcal{R}_{t},\mathbf{0 }\right]+\left[\mathcal{R}_{t},\mathcal{R}^{\prime}_{t}\right], \tag{5}\]

where \(\mathbf{\Theta}_{t}(\mathcal{X}^{\prime},\mathcal{X})\triangleq\nabla_{\mathbf{W} }\mathcal{F}^{\top}_{t}\nabla_{\mathbf{W}}\mathcal{F}_{t}\in\mathbb{R}^{(n-n )\times n_{l}}\) is the NTK matrix between training and testing sets, which quantifies similarity between seen and unseen instances based on how differently their outputs change by an infinitesimal perturbation of weights. The whole \(n\times n\) NTK matrix will be henceforth abbreviated as \(\mathbf{\Theta}_{t}\). This equation implies that the GD-based optimization of arbitrary parameterized models can be viewed as an information propagation process where residuals propagate unidirectionally from training to arbitrary unseen samples based on a kernel similarity measure.

**Residual Propagation.** Drawing upon an analogy between the information propagation process in (5) induced by optimization, and message passing schemes between instances commonly seen in graph learning, we propose to replace the dense NTK matrix \(\mathbf{\Theta}_{t}\) with a sparse matrix (e.g. high-order graph adjacency matrix \(\mathbf{A}^{K}\) ) as a similarity measure, which allows us to implement (5) as a new semi-supervised algorithm called _Residual Propagation (RP)_ with update equation:

\[\left[\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\right]=-\eta\mathbf{A}^{K} [\mathcal{R}_{t},\mathbf{0}]+[\mathcal{R}_{t},\mathcal{R}^{\prime}_{t}],\ \text{ where }\mathcal{R}_{0}=\mathcal{Y}, \mathcal{R}^{\prime}_{0}=\mathbf{0}. \tag{6}\]

As initialization, the model predictions \(\mathcal{F}_{0}\) and \(\mathcal{F}^{\prime}_{0}\) are defined as \(\mathbf{0}\), and unknown testing labels are defined as \(\mathcal{Y}^{\prime}=\mathbf{0}\) such that negative \(\mathcal{R}^{\prime}_{t}\) equals model predictions. Compared with (5) merely as an expression, RP is a practical algorithm without trainable weights or relying on input features, and is efficient with linear time complexity owing to sparse matrix multiplication. For further intuition, consider an arbitrary unseen data point \(\mathbf{x}^{\prime}\), whose update equation in RP is

\[\boxed{\mathbf{\frac{\mathbf{\mathsf{P}_{t}}}{\mathbf{\mathsf{f}_{t+1}(\mathbf{x}^{\prime})}} }=f_{t}(\mathbf{x}^{\prime})+\eta\sum_{\mathbf{x}\in\mathcal{X}}\mathbf{A}^{K}_{\mathbf{x} ^{\prime}\mathbf{x}}\ \left(\begin{array}{c}\frac{1}{\left(\mathbf{y}(\mathbf{x})\right)}-f_{t}(\mathbf{x}) \end{array}\right), \tag{7}\]

where \(y(\mathbf{x})\) is the ground-truth label for \(\mathbf{x}\), and \(\mathbf{A}^{K}_{\mathbf{x}^{\prime}\mathbf{x}}\) is an element of \(\mathbf{A}^{K}\) corresponding to the sample pair \((\mathbf{x},\mathbf{x}^{\prime})\). For an unseen instance \(\mathbf{x}^{\prime}\) that is more'similar' or 'close' to \(\mathbf{x}\), more ground-truth label information \(y(\mathbf{x})\) will then propagate to the testing instance, and vice versa, to shape the prediction for the unseen one (illustrated in Fig. 1(a)). For \(f_{t}(\mathbf{x})\neq 0\), the ground-truth label is adjusted by subtracting current model prediction, i.e. \(y(\mathbf{x})-f_{t}(\mathbf{x})\), enabling the propagation process to diminish progressively as the predictions become more accurate.

### Theoretical Analysis and Empirical Evaluation

Intriguingly, we found the RP algorithm is connected with various classic methods including LP and kernel regression, though they emerge from very different contexts. (Analysis here for RP is also useful for understanding GNNs as will be discussed in the next section.)

**Proposition 1** (Connection with Label Propagation).: _The first step of RP in (6) yields identical classification results as LP in (3) (with \(\alpha=1\) and \(k=K\)):_

\[\text{(First Step of RP): }[\mathcal{F}_{1},\mathcal{F}^{\prime}_{1}]=\eta\mathbf{A}^{K}[ \mathcal{Y},\mathbf{0}],\ \ \text{(Label Propagation): }\ \ \text{LP}(\mathcal{Y};K,1)=\mathbf{A}^{K}[\mathcal{Y},\mathbf{0}]. \tag{8}\]

_Remark_.: Besides the first step, each of subsequent step of RP can also be viewed as LP on adjusted ground-truth labels, i.e. \(\mathcal{Y}-\mathcal{F}_{t}=\mathcal{R}_{t}\). This connection offers the potential for improvement as we iterate beyond the initial LP solution, and also demonstrates the flexibility of RP by replacing
with a more general function \(\operatorname{LP}^{*}(\cdot):\mathbb{R}^{n_{l}}\to\mathbb{R}^{n}\) that takes as input ground-truth labels and outputs predictions, leading to generalized RP with the following update equation:

\[\left[\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\right]=-\eta \operatorname{LP}^{*}(\mathcal{R}_{t})+\left[\mathcal{R}_{t},\mathcal{R}^{ \prime}_{t}\right]. \tag{9}\]

**Theorem 2** (Convergence and Optimization Guarantee).: _For RP in (6) and sufficiently small step size \(\eta<2\sigma_{max}^{-1}[\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}]\), where \(\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}\) is a submatrix of \(\mathbf{A}^{K}\), and \(\sigma_{max}[\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}]\) is its largest eigenvalue, \(\mathcal{R}_{t}\) and \(\mathcal{R}^{\prime}_{t}\) converge as \(t\to\infty\) for positive definite \(\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}\) or positive semi-definite \(\mathbf{A}^{K}\)._

**Corollary 3** (Connection with Kernel Regression).: _Upon convergence for positive definite \(\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}\), the model predictions are equivalent to the kernel regression solution w.r.t. kernel \(\kappa(\mathbf{x},\mathbf{x}^{\prime})\triangleq\mathbf{A}_{\mathbf{x}\mathbf{x}^{\prime}}^{K}\)_

\[\mathcal{F}_{\infty}=\mathcal{Y}\quad\text{(Perfect fit of training set)},\quad\mathcal{F}_{\infty}^{\prime}=\mathbf{A}_{\mathcal{X}\mathcal{X}}^{K}( \mathbf{A}_{\mathcal{X}\mathcal{X}}^{K})^{-1}\mathcal{Y}. \tag{10}\]

_Remark_.: There exists one exception though where the convergence is not guaranteed, namely when \(\mathbf{A}\) has negative eigenvalues and \(K\) is odd. While this can be easily avoided by adding diagonal components (e.g. PageRank matrix \(\alpha\mathbf{A}+(1-\alpha)\mathbf{I}_{n}\)) for enforcing positive semi-definiteness, empirically we do not need the algorithm to converge by stopping at the step with peak validation performance (just like standard training procedure for neural networks), and the algorithm can still achieve satisfactory generalization performance (often better than the version of RP that can converge).

**Empirical Evaluation.** We compare the basic version of RP from (6) with some representative GNN architectures (LinearGNN Wu et al. (2019) and GCN Kipf and Welling (2017)) on challenging OGB (Hu et al., 2020) node classification datasets Arxiv, Proteins, Products with up to millions of nodes and edges. Details of experimental setting are deferred to Appendix D.1. The codes will be available at [https://github.com/chr26195/ResidualPropagation](https://github.com/chr26195/ResidualPropagation).

\(\circ\)_Effectiveness:_ In Table 1, the proposed RP demonstrates decent performance, often surpassing GNNs even when no node or edge feature information is utilized. Notably, this performance even approaches some other more advanced models such as node-level graph Transformers (Wu et al., 2022). As depicted in Fig. 2, RP achieves the same performance as LP using one step, and quickly increases until reaching its peak performance, which surpasses the GNN. Specifically, on Proteins where the graph contains relatively richer structural information, a single step of RP exceeds a well-trained deep GNN, while in Products, \(4\) steps of RP exceeds the GNN.

\(\circ\)_Efficiency and Scalability:_ In addition to its effectiveness, RP offers several advantages. It does not require learnable parameters and boosts speed, with each step being more than 10 times faster than the GNN, in part because it does not need a backward pass typically required in conventional

\begin{table}
\begin{tabular}{c|c|c c|c c|c c|c|c} \hline \hline
**Model** & **Feat.** & \multicolumn{2}{c|}{Arxiv} & \multicolumn{2}{c|}{Proteins} & \multicolumn{2}{c|}{Products} & \multicolumn{2}{c|}{**Avg.**} & \multicolumn{2}{c}{**Param.**} \\  & & Validation & Test & Validation & Test & Validation & Test & \multicolumn{2}{c|}{**Test**} & \\ \hline MLP & \(\mathcal{X}\) & \(57.65\pm 0.12\) & \(55.50\pm 0.23\) & \(77.06\pm 0.14\) & \(72.04\pm 0.48\) & \(75.54\pm 0.14\) & \(61.06\pm 0.08\) & \(66.48\) & \(O(dm^{2})\) \\ LinearGNN & \(\mathcal{X}\), \(\mathbf{A}\) & \(70.67\pm 0.02\) & \(69.39\pm 0.11\) & \(66.11\pm 0.87\) & \(62.99\pm 0.11\) & \(88.97\pm 0.01\) & \(74.21\pm 0.04\) & \(72.04\) & \(O(dc)\) \\ GNN & \(\mathcal{X}\), \(\mathbf{A}\) & \(\mathbf{730.0\pm 0.17}\) & \(71.74\pm 0.29\) & \(79.21\pm 0.18\) & \(72.51\pm 0.03\) & \(92.00\pm 0.00\) & \(75.64\pm 0.021\) & \(72.35\) & \(O(dm^{2})\) \\ LP & \(\mathbf{A}\) & \(70.14\pm 0.08\) & \(68.32\pm 0.00\) & \(83.02\pm 0.00\) & \(76.23\pm 0.00\) & \(90.91\pm 0.00\) & \(74.34\pm 0.00\) & \(76.91\) & \(0\) \\ \hline RP (ours) & A & \(71.37\pm 0.00\) & \(70.06\pm 0.00\) & \(85.19\pm 0.00\) & \(78.17\pm 0.00\) & \(91.31\pm 0.00\) & \(78.25\pm 0.00\) & \(79.06\) & \(0\) \\ Speedup / step & & \(\times\) 14.48 & \(\times\) 14.00 & \(\times\) 12.46 & & & & \\ Time to Acc. & & \(\times\) 0.01461 & & \(\times\) 0.0008 & \(\times\) 0.00427 & & & \\ Memory & & \(\times\) 0.094 & & \(\times\) 0.363 & & \(\times\) 0.151 & & & \\ \hline \hline \end{tabular}
\end{table}
Table 1: Empirical evaluation of RP on OGB datasets. Accuracy is reported for Arxiv and Products, and ROC-AUC for Proteins. Last three rows compare RP against full-batch GNN.

Figure 2: Learning curves of RP and comparison with the performance of LP (\(\alpha=1\)), linear GNN and deep GNN. Transition from yellow to purple denotes RP with decreasing step size \(\eta\).
learning frameworks. Furthermore, RP achieves its peak performance within a mere \(20\) steps, and thus overall takes orders of magnitudes less time (\(<0.01\)) to attain GNN-level performance. In terms of memory cost, RP inherits the scalability of LP and only requires storage space for model predictions; in practice, we often found the memory bottleneck lies with data preprocessing \(O(|\mathcal{E}|)\) rather than the algorithm itself \(O(n)\).

In AppendicesE.1 andE.2 respectively, we discuss extensions of RP that combine with kernels (e.g. Gaussian kernel) to accommodate node features, and test additional \(7\) datasets (Cora, Citeseer, Pubmed, Computer, Photo, CS, Physics) where RP can still outperform popular GNNs.

## 4 Learning Dynamics of Graph Neural Networks

We have introduced a simple propagation method that competes with deep GNN models, which raises the question whether or not GNNs are also (secretly) taking advantage of similar residual propagating effects with a graph inductive bias in their optimization process. To examine this possibility, we formally study the learning dynamics of GNNs in function space during GD-based training.

### General Form by Node-Level GNTK

Similar to the learning dynamics of general parameterized models in (5), the learning dynamics of GNNs in node-level tasks is characterized by their NTK defined as follows:

**Definition 4** (Node-Level Graph Neural Tangent Kernel).: _For a \(\ell\)-layer GNN in node-level tasks defined in Sec. 2, the NTK is defined as_

\[\mathbf{\Theta}_{t}^{(\ell)}(\mathbf{x},\mathbf{x}^{\prime};\mathbf{A})=\nabla_{\mathbf{W }}f(\mathbf{x};\mathbf{A})^{\top}\nabla_{\mathbf{W}}f(\mathbf{x}^{\prime};\mathbf{A}), \tag{11}\]

_which we refer to as Node-Level Graph Neural Tangent Kernel (GNTK) to differentiate it with the graph-level GNTK for graph-level tasks proposed in Du et al. (2019), or simply NTK of GNNs._

How GNNs evolve in function space during training also follows the residual propagation process:

\[\begin{bmatrix}\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\end{bmatrix}=- \eta\mathbf{\Theta}_{t}^{(\ell)}[\mathcal{R}_{t},\mathbf{0}]+[\mathcal{R}_{t}, \mathcal{R}^{\prime}_{t}]\,,\quad\mathbf{\Theta}_{t}^{(\ell)}\triangleq[\mathbf{ \Theta}_{t}^{(\ell)}(\mathbf{x}_{i},\mathbf{x}_{j};\mathbf{A})]_{i,j\in[n]\times[n]}, \tag{12}\]

where the initial values of residuals depend on initial weights. A precise characterization of (12) requires mathematically computing the node-level GNTK, which is prohibitively challenging for finitely-wide GNNs. Fortunately, we can derive its explicit formula in overparameterized regimes where the model width \(m\) tends to infinity, following that of fully-connected neural networks (Jacot et al., 2018; Lee et al., 2019). Given that the limiting NTK in this regime behaves as a constant kernel function (Yang and Littwin, 2021), we change the notation by omitting the subscript \(t\).

Concretely, we present the computation of the node-level GNTK in a recurrent layer-wise manner where each layer is decomposed into two steps: _Feature Transformation_ (corresp. \(\mathbf{Z}\leftarrow\sigma(\mathbf{Z}\mathbf{W})\)) and _Feature Propagation_ (corresp. \(\mathbf{Z}\leftarrow\mathbf{A}\mathbf{Z}\)). As elements in the computation, we denote a node-level GNTK for GNN without feature propagation at \(\ell\)-th layer as \(\mathbf{\Theta}^{(\ell)}\), the covariance matrix of the \(\ell\)-th layer's outputs with (without) feature propagation as \(\mathbf{\Sigma}^{(\ell)}\) (\(\mathbf{\bar{\Sigma}}^{(\ell)}\)), and the covariance matrix of the derivative to the \(\ell\)-th layer as \(\mathbf{\bar{\Sigma}}^{(\ell)}\). Then the feature propagation and transformation steps in each layer respectively correspond to (we concisely show the key steps here and defer the complete formula to AppendixB.3):

\[\begin{array}{ll}\text{(Feature Transformation)}&\text{(Feature Propagation)}\\ &\\ \mathbf{\bar{\Theta}}^{(\ell)}=\mathbf{\Theta}^{(\ell-1)}\odot\mathbf{\bar{\Sigma}}^{(\ell )}+\mathbf{\bar{\Sigma}}^{(\ell)},&\\ &\\ \mathbf{\Theta}^{(\ell)}=\mathbf{A}\ \mathbf{\bar{\Theta}}^{(\ell)}\mathbf{A}.\end{array} \tag{13}\]

**Implications.** Compared with the NTK computation for a fully-connected neural network (Jacot et al., 2018), the node-level GNTK has an equivalent feature transformation step, while its uniqueness stems from the feature propagation step, whereby the adjacency matrix \(\mathbf{A}\) integrates into the kernel similarity measure. Consequently, this kernel function naturally accommodates a graph inductive bias, and thus the residual propagation process of GNNs also tends to follow the trajectory regulated by the graph, similar to the behavior of the RP algorithm. While such an analysis offers valuable mathematical intuitions, we next present how the learning dynamics of certain GNNs in function space can be exactly formulated under the framework of generalized RP.
### Case Studies: from Two Layers to Arbitrary Depth

Given that our primary focus centers on the role of graphs in GNNs (as without them, GNNs are largely equivalent to MLPs), we exclude external node features and instead define inputs as either: 1) an identity matrix \(\bar{\mathcal{X}}\triangleq\mathbf{I}_{n}\) that assigns each node a one-hot vector as indication of its unique identity (as sometimes assumed in practice (Kipf and Welling, 2017; Zhu et al., 2021)), which can be viewed as learning a unique embedding for each node by treating the first-layer weights as an embedding table; 2) fixed node embeddings from graph spectral decomposition \(\bar{\mathcal{X}}\triangleq\operatorname*{arg\,min}_{\mathbf{B}}\|\mathbf{A}- \mathbf{B}\mathbf{B}^{\top}\|_{F}^{2}\), which aligns with various network embedding approaches based on definitions of \(\mathbf{A}\)(Qiu et al., 2018).

**Two-Layer GNN.** Building on the setup from prior work analyzing infinitely-wide two-layer neural networks (Arora et al., 2019), we consider two-layer GNNs:

**Theorem 5** (Two-Layer GNN Dynamics).: _For an infinitely-wide two-layer GNN defined as \([\mathcal{F},\mathcal{F}^{\prime}]=\mathbf{A}\sigma(\mathbf{A}\bar{\mathcal{ X}}\mathbf{W}^{(1)})\mathbf{W}^{(2)}/\sqrt{m}\) with ReLU activation, \(\bar{\mathcal{X}}=\mathbf{I}_{n}\) and standard NTK parameterization, its learning dynamics by optimizing \(\mathbf{W}^{(1)}\) can be written as a generalized RP process_

\[\left[\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\right]=-\eta\mathbf{A}( \mathbf{A}^{2}\odot\mathbf{S})\mathbf{A}[\mathcal{R}_{t},\mathbf{0}]+[ \mathcal{R}_{t},\mathcal{R}^{\prime}_{t}],\ \ \mathbf{S}_{ij}=\left(\pi-\arccos(\frac{ \mathbf{A}^{\top}\mathbf{A}_{j}}{\|\mathbf{A}_{i}\|\mathbf{A}_{j}\|})\right) /2\pi. \tag{14}\]

_The matrix \(\mathbf{S}\) reweights each element in \(\mathbf{A}^{2}\) by the similarity of neighborhood distributions of two nodes. For \(\bar{\mathcal{X}}=\operatorname*{arg\,min}_{\mathbf{B}}\|\mathbf{A}-\mathbf{B }\mathbf{B}^{\top}\|_{F}^{2}\), the propagation matrix is replaced by \(\mathbf{A}(\mathbf{A}^{3}\odot\tilde{\mathbf{S}})\mathbf{A}\) where \(\tilde{\mathbf{S}}\) is another reweighting matrix (details given in Appendix B.4)._

**Arbitrarily Deep GNN.** Pushing further, we can also characterize the evolution of arbitrarily deep GNNs where feature propagation is applied at the last layer (e.g. Klicpera et al. (2019); Liu et al. (2020); Spinelli et al. (2020); Chien et al. (2021), called decoupled GNNs in some other literature):

**Theorem 6** (Deep and Wide GNN Dynamics).: _For arbitrarily deep and infinitely-wide GNNs with feature propagation deferred to the last layer, i.e. \([\mathcal{F},\mathcal{F}^{\prime}]=\mathbf{A}^{\ell}\operatorname{MLP}(\bar{ \mathcal{X}})\) with \(\bar{\mathcal{X}}=\mathbf{I}_{n}\), the learning dynamics that result from optimizing MLP weights can be written as the generalized RP process_

\[\left[\mathcal{R}_{t+1},\mathcal{R}^{\prime}_{t+1}\right]=-\eta\mathbf{A}^{ \ell}(\mathbf{I}_{n}+c\mathbf{1}\mathbf{1}^{\top})\mathbf{A}^{\ell}[\mathcal{R }_{t},\mathbf{0}]\ +\ \left[\mathcal{R}_{t},\mathcal{R}^{\prime}_{t}\right], \tag{15}\]

_where \(c\geq 0\) is a constant determined by the model depth, and \(\mathbf{1}\) is an all-\(1\) column vector._

As a special case of the above, when the backbone MLP only has one layer, the model degrades to a linear GNN (e.g. SGC Wu et al. (2019)), and the constant correspondingly is \(c=0\), giving us the basic version of RP in (6) or its approximation for odd \(K\):

**Corollary 7** (Linear GNN Dynamics).: _The learning dynamics of the linear GNN \([\mathcal{F},\mathcal{F}^{\prime}]=\mathbf{A}^{\ell}\bar{\mathcal{X}}\mathbf{W}\) is identical to the basic version of RP in (6) with \(K=2\ell\) for input features \(\bar{\mathcal{X}}=\mathbf{I}_{n}\), and \(K=2\ell+1\) for \(\bar{\mathcal{X}}=\operatorname*{arg\,min}_{\mathbf{B}}\|\mathbf{A}-\mathbf{B }\mathbf{B}^{\top}\|_{F}^{2}\) and positive semi-definite \(\mathbf{A}\)._

_Remark_.: Despite this equivalence, it is important to note that this specific linear GNN (on our given input features) is not as lightweight as it may appear. First, its parameter number scales with the size of dataset, often reaching orders of magnitude larger than deep GNN models (e.g., \(115,104,363\) for the linear model v.s. \(103,727\) for GCN on Products). Additionally, the full-rank spectral decomposition of \(\mathbf{A}\) is computationally very expensive (\(O(n^{3})\)) for large graphs. In contrast, the basic version RP in (6) efficiently yields identical results to this heavily parameterized GNN without actually training parameters or computing expensive matrix decompositions. Furthermore, the generalized RP from (9) can implement algorithms that are infeasible within a conventional deep learning framework, e.g. training infinitely-wide GNNs.

## 5 Rethinking the Success and Pathology of GNNs

In the last section, we have given theoretical analysis with concrete examples to illustrate that NTKs of GNNs naturally tend to align with the graph. This finding, coupled with the consistently strong performance of the RP algorithm -- an example of perfect graph-kernel alignment -- prompts the question of whether this alignment between NTK and the graph plays a crucial role in understanding the success of GNNs. Next, we offer interpretable explanations of "when and why GNNs successfully generalize" and their pathological training behavior on heterophilic graphs (Sec. 5.1). We also study the time evolution of real-world GNN NTKs to empirically verify our theoretical results (Sec. 5.2).
### When and Why GNNs Successfully Generalize?

Our previous discussions have revolved around two matrices, namely the graph adjacency matrix \(\mathbf{A}\) and the NTK matrix \(\mathbf{\Theta}_{t}\).1 To complete the theoretical picture, we introduce another matrix called the ideal _or optimal kernel matrix_(Cristianini et al., 2001), defined as \(\mathbf{\Theta}^{*}\triangleq\tilde{\mathcal{Y}}\tilde{\mathcal{Y}}^{\top}\in \mathbb{R}^{n\times n}\) to indicate whether two instances have the same label, and a metric to quantify alignment of (kernel) matrices:

Footnote 1: We here refer to \(\mathbf{A}\) as a class of similarity matrices based on original \(\mathbf{A}\) in a general sense, such as \(\mathbf{A}^{K}\) etc.

**Definition 8** (Alignment, Cristianini et al. (2001)).: _Given two (kernel) matrices \(\mathbf{K}_{1}\) and \(\mathbf{K}_{2}\), their alignment is defined as \(A\left(\mathbf{K}_{1},\mathbf{K}_{2}\right)\triangleq\left(\mathbf{K}_{1}, \mathbf{K}_{2}\right)_{F}/(\|\mathbf{K}_{1}\|_{F}\|\mathbf{K}_{2}\|_{F})\in[0,1]\). This is a generalization of cosine similarity from vectors to matrices, \(\arccos\) of which satisfies the triangle inequality._

In the following, we will first take a separate look at the effects of alignment for each pair of the three matrices of interest (cf. Fig. 1(a)), and then discuss their collective implications.

\(\circ\)_Homophily Level:_\(A(\mathbf{A},\mathbf{\Theta}^{*})\). The alignment between \(\mathbf{A}\) and \(\mathbf{\Theta}^{*}\) quantifies the _homophily level_ of graph structure, i.e. whether two connected nodes indeed have the same label, and is determined the dataset. While many empirical results (e.g. Zhu et al. (2020)) suggest high homophily level is important for the performance of GNNs, deeper theoretical understandings are mostly lacking.

\(\circ\)_Kernel-Target Alignment:_\(A(\mathbf{\Theta}_{t},\mathbf{\Theta}^{*})\). The alignment between kernel matrix and optimal \(\mathbf{\Theta}^{*}\) has been widely studied in the research on kernel learning (Cristianini et al., 2001; Kwok and Tsang, 2003; Lanckriet et al., 2004; Gonen and Alpaydin, 2011). Better kernel-target alignment has been recognized as a critical factor that leads to favorable generalization. For intuition, one can quickly verify that substituting \(\mathbf{\Theta}^{*}\) to the residual dynamics in (5) leads to perfect generalization performance (since ground-truth labels only propagate to unseen instances with the same label).

\(\circ\)_Kernel-Graph Alignment:_\(A(\mathbf{\Theta}_{t},\mathbf{A})\). The alignment between NTK and graph is a novel notion in our work, as prior sections have shown that GNN NTK matrices naturally tend to align with \(\mathbf{A}\). E.g., the RP algorithm (and variants thereof) serve as an extreme case with two identical matrices.

**Implications.** We consider two cases. For _homophilic_ graphs where \(A(\mathbf{A},\mathbf{\Theta}^{*})\uparrow\) is naturally large, better kernel-graph alignment \(A(\mathbf{\Theta}_{t},\mathbf{A})\uparrow\) consequently leads to better kernel-target alignment \(A(\mathbf{\Theta}_{t},\mathbf{\Theta}^{*})\uparrow\). In other words, the NTK of GNNs naturally approaches the optimum as the graph structure possesses homophily property, and leveraging it in the optimization process (12) encourages training residuals to flow to unseen samples with the same label and thus better generalization; In contrast, for _heterophilic_ graphs where \(A(\mathbf{A},\mathbf{\Theta}^{*})\) is small, better kernel-graph alignment will hinder kernel-target alignment, explaining the pathological learning behavior of GNNs when dealing with heterophilic graphs in an interpretable manner.

**Theoretical Support 1.** To support this interpretation, we examine the generalization behavior of infinitely-wide neural networks in the extreme case where its NTK matrix is assumed to be governed by the graph, say, \(\lim_{k\to\infty}\sum_{i=0}^{k}(\alpha\mathbf{A})^{i}\) with \(\alpha\in(0,1)\) as adopted by the converged LP algorithm in (3). With a common assumption that training instances are drawn i.i.d. from a distribution \(\mathcal{P}\), and based on the Rademacher complexity generalization bound for kernel regression (Bartlett and Mendelson, 2002; Arora et al., 2019; Du et al., 2019), we have a label-dependent high-probability (at least \(1-\delta\)) upper bound on population risk (derivation in Appendix B.7):

\[\mathbb{E}_{(\mathbf{x},y)\sim\mathcal{P}}\left[l\left(f(\mathbf{x}),y\right)\right]= O\left(\sqrt{\frac{n_{l}-cA(\mathbf{A},\mathbf{\Theta}^{*})}{n_{l}}}+\sqrt{ \frac{\log(1/\delta)}{n_{l}}}\right), \tag{16}\]

where \(c=\alpha\|\mathbf{\Theta}^{*}\|_{F}\|\mathbf{A}\|_{F}\) is a constant, \(A(\mathbf{A},\mathbf{\Theta}^{*})\) is the homophily level for the training set (with slight abuse of notation). This bound can also be viewed as a theoretical guarantee for the converged (generalized) RP algorithm and clearly demonstrates that: _for\(A(\mathbf{\Theta}_{t},\mathbf{A})\) fixed as \(1\), higher level of graph homophily \(A(\mathbf{A},\mathbf{\Theta}^{*})\) plays a dominant role in better generalization_. The above analysis can also be potentially extended to infinitely-wide GNNs discussed in Sec. 4 with further assumptions on the condition number of its NTK matrix, which will lead to similar (but less straightforward) results.

**Theoretical Support 2.** Pushing further, in a complementary setting to the above, where the objective is to find the optimal a priori kernel for directly minimizing the population risk (but without access to any input features or ground-truth labels), we demonstrate that when homophily assumptions are
imposed on the graph, infinitely-wide GNNs and the RP algorithm will yield the optimal kernel regression predictor with the provably best generalization performance:

**Theorem 9** (Bayesian Optimality of GNN).: _Assume that the underlying data generation distribution \(\mathcal{P}\) is such that the probability of a pair of instances having the same label \(P(y(\mathbf{x}_{i})=y(\mathbf{x}_{j}))\) is proportional to \(\mathbf{A}_{ij}-1/2\). Then the optimal kernel regression predictor that minimizes the population risk with squared loss has kernel matrix \(\mathbf{A}\)._

_Remark_.: The matrix \(\mathbf{A}\) from the above result could vary depending on different assumptions on \(P(y(\mathbf{x}_{i})=y(\mathbf{x}_{j}))\), such as \(\mathbf{A}^{K}\) for the basic version of RP, or \(\mathbf{A}\left(\mathbf{A}^{2}\odot\mathbf{S}\right)\mathbf{A}\) for infinitely-wide two-layer GNNs in Theorem 5. And ideally, if one has privileged access to labels of all instances, this optimal matrix is exactly the optimal (posterior) kernel \(\mathbf{\Theta}^{*}\) discussed above and perfect generalization will be achieved. This result further verifies our interpretation of generalization by showing that _for \(A(\mathbf{A},\mathbf{\Theta}^{*})\) fixed to be large, better \(A(\Theta_{\ell},\mathbf{A})\) (e.g. GNN and RP) leads to favorable generalization.

### Empirical Verification and Further Insights on Heterophily

We next empirically verify results from Sec. 4 and 5.1 by studying the time evolution of real-world GNN NTKs. Figure 3 plots the results reflected by alignment between \(\mathbf{\Theta}_{t}\), \(\mathbf{A}\) and \(\mathbf{\Theta}^{*}\).

**Synthetic Dataset.** To demonstrate the effects of different homophily levels, we use a stochastic block model (Holland et al., 1983) to generate a synthetic dataset with informative node features. As shown in Fig. 3(a), the kernel-graph alignment \(A(\Theta_{t},\mathbf{A})\) for GNNs stays at a high level regardless of different graph structures, as a natural result of the network architecture (Sec. 4). Consequently, as we gradually alter edge probabilities to decrease the homophily level \(A\left(\mathbf{A},\mathbf{\Theta}^{*}\right)\), the kernel-target alignment \(A\left(\Theta_{t},\mathbf{\Theta}^{*}\right)\) also decreases and the testing accuracy drops dramatically (Sec. 5.1).

**Real-World Datasets.** On real-world homophilic and heterophilic datasets, we progressively coarsen the graph until the GNN degrades to an MLP, allowing us to analyze the effects of feature propagation in the model. For now, let us first focus on comparing red and blue lines in Fig. 3(b,c). We found the kernel-graph alignment \(A(\Theta_{t},\mathbf{A})\) overall decreases with less graph structure, again verifying results in Sec. 4. However, the impact of the graph structure depends on different homophily levels \(A\left(\mathbf{A},\mathbf{\Theta}^{*}\right)\): on the homophilic dataset, more graph structure and better \(A(\mathbf{\Theta}_{t},\mathbf{A})\) optimize the NTK matrix as reflected by better kernel-target alignment \(A\left(\mathbf{\Theta}_{t},\mathbf{\Theta}^{*}\right)\), but worsens it on the heterophilic one. This results in distinct generalization behavior of GNNs reflected by the testing accuracy.

**Can and How (G)NNs Handle Heterophily?** Recently, there has been a growing discussion on whether standard GNNs are capable of handling heterophily, with empirical results pointing to diverging conclusions (Zhu et al., 2020; Ma et al., 2022; Luan et al., 2022; Platonov et al., 2023). From the learning dynamics perspective, we add new insights to this debate: \(\circ\)**1)** While we have shown that infinitely-wide GNNs are sub-optimal on heterophilic graphs, MLPs in contrast have

Figure 3: Evolution of NTK matrix \(\mathbf{\Theta}_{t}\) of finitely-wide GCN during training, reflected by matrix alignment. **(a)** Synthetic dataset generated by a stochastic block model, where the homophily level gradually decreases by altering edge probabilities, i.e. homophilic \(\rightarrow\) heterophilic; **(b & c)** Real-world homophilic (Cora) and heterophilic (Texas) datasets, where the graph is gradually coarsened until there is no edge left when evaluating \(\mathbf{\Theta}_{t}\), i.e. more graph \(\rightarrow\) less graph. (Details in Appendix D.2)
no guarantee of kernel-target alignment (since both \(A(\mathbf{A},\mathbf{\Theta}^{*})\) and \(A(\mathbf{\Theta}_{t},\mathbf{A})\) are not well-aligned), indicating that _they could either generalize better or worse without clearcut (dis)advantages, as opposed to the case on homophilic graphs where GNN is provably at an advantage_, which explains the diverse empirical results when comparing GNNs with MLPs in prior work. \(\circ\)**2)** As we now turn to the NTK's time evolution in Fig. 3, an additional phenomenon we found across all datasets is the overall increase of kernel-target alignment \(A\left(\mathbf{\Theta}_{t},\mathbf{\Theta}^{*}\right)\) during the training process. This indicates that the training process enables real-world GNNs to adjust their NTK's feature space such that the kernel matrix learns towards an ultimately homophilic graph structure (i.e. \(\mathbf{\Theta}^{*}\)) to adapt to heterophilic datasets (and consequently different evolutionary trends of \(A\left(\mathbf{\Theta}_{t},\mathbf{A}\right)\) for homophilic and heterophilic datasets). Such a phenomenon has also been found for other models in vision tasks (Baratin et al., 2021). Additionally, since our analysis also applies to other forms of \(\mathbf{A}\) for feature propagation in GNNs, it could also potentially explain how some specialized models with different definitions of \(\mathbf{A}\) (e.g. allowing signed propagation) can mitigate the heterophily issue.

In Appendix E.3, we consider additional \(5\) heterophilic datasets and compare RP against recently proposed GNNs specially designed for handing heterophily. We found this simple algorithm still achieves decent performance despite it is theoretically sub-optimal.

## 6 More Discussions and Conclusion

**Abridged Related Work.** Most existing work on theoretical aspects of GNNs focuses on representation and generalization (see Jegelka (2022) and references therein) while their optimization properties remains under-explored (Zhang et al., 2020; Xu et al., 2021; Yadati, 2022). Specifically, existing work in representation (or expressiveness) does not provide answers to what exactly GNN functions are found during the optimization process. For generalization, prior work is insufficient to explain the effects of training, which is widely-recognized as a crucial ingredient, and does not connect to heterophily, a relevant aspect in practice. (See unabridged related work in Appendix A.)

**Applicability.** Our insights apply to both transductive and inductive settings, other loss functions, and multi-dimensional outputs (see Appendix C). The analytical framework could also be extended to other tasks (which are left as future work): for graph classification or regression, the definition of GNTK in (12) should be modified to the graph-level one (Du et al., 2019); for link prediction and self-supervised learning, the loss function in the derivation of (5) should be adjusted accordingly. A similar message passing process in function space would still hold in these settings. Moreover, the proposed RP is compatible with regression tasks due to minimization of the squared loss.

**Broader Impacts.** Theoretically, our framework of interpreting GD dynamics in function space as a label propagation process across seen and unseen samples could be adapted for dissecting other models and tasks, and the analysis in Sec. 5 exemplifies using kernel learning to justify model architectures. Practically, the proposed RP unlocks the possibility of leveraging a sparse structure in complicated training dynamics, acting as efficient forward-only solutions to implement dynamics in function space which are previously impractical or extremely expansive (Arora et al., 2020). More broadly, it also serves as a bridge of classic graph-based algorithms (Chapelle et al., 2009) and modern deep learning, and echoes recent efforts to develop forward-only learning frameworks (Hinton, 2022).

**Conclusion.** This work represents a preliminary endeavor in formally studying the learning dynamics of GNNs in function space. We reveal the alignment between external graph structure and NTK matrix that controls the evolution of models, which enables deeper understanding of the generalization capability and pathological behavior of GNNs by connecting to the study of heterophily, and leads to a minimalist algorithm with superior efficiency, decent performance and theoretical interpretability.

## Ethics Statement

In this research, we investigate the learning dynamics of GNNs in function space and offer insights to their optimization and generalization behaviors. Our insights can deepen understandings for graph-based deep learning, and do not appear to carry immediate harmful consequences to the best of our knowledge.
## Reproducibility Statement

The complete proofs of our theoretical results and their assumptions are provided in Appendix B. Extensive implementation details of the experiments from Section 3.2 and Section 5.2 regarding baselines, hyperparameters and datasets are provided in Appendix D.1 and D.2. Detailed descriptions for the implementation of different versions of RP algorithms are given in Appendix D.3.
