Neural Tangent Kernels Motivate Graph Neural Networks with
  Cross-Covariance Graphs
Neural tangent kernels (NTKs) provide a theoretical regime to analyze the
learning and generalization behavior of over-parametrized neural networks. For
a supervised learning task, the association between the eigenvectors of the NTK
kernel and given data (a concept referred to as alignment in this paper) can
govern the rate of convergence of gradient descent, as well as generalization
to unseen data. Building upon this concept, we investigate NTKs and alignment
in the context of graph neural networks (GNNs), where our analysis reveals that
optimizing alignment translates to optimizing the graph representation or the
graph shift operator in a GNN. Our results further establish the theoretical
guarantees on the optimality of the alignment for a two-layer GNN and these
guarantees are characterized by the graph shift operator being a function of
the cross-covariance between the input and the output data. The theoretical
insights drawn from the analysis of NTKs are validated by our experiments
focused on a multi-variate time series prediction task for a publicly available
dataset. Specifically, they demonstrate that GNNs with cross-covariance as the
graph shift operator indeed outperform those that operate on the covariance
matrix from only the input data.
Shervin Khalafi, Saurabh Sihag, Alejandro Ribeiro
# Neural Tangent Kernels Motivate Graph Neural Networks with Cross-Covariance Graphs

###### Abstract

Neural tangent kernels (NTKs) provide a theoretical regime to analyze the learning and generalization behavior of over-parametrized neural networks. For a supervised learning task, the association between the eigenvectors of the NTK kernel and given data (a concept referred to as alignment in this paper) can govern the rate of convergence of gradient descent, as well as generalization to unseen data. Building upon this concept, we investigate NTKs and alignment in the context of graph neural networks (GNNs), where our analysis reveals that optimizing alignment translates to optimizing the graph representation or the graph shift operator in a GNN. Our results further establish the theoretical guarantees on the optimality of the alignment for a two-layer GNN and these guarantees are characterized by the graph shift operator being a function of the cross-covariance between the input and the output data. The theoretical insights drawn from the analysis of NTKs are validated by our experiments focused on a multi-variate time series prediction task for a publicly available dataset. Specifically, they demonstrate that GNNs with cross-covariance as the graph shift operator indeed outperform those that operate on the covariance matrix from only the input data.

## 1 Introduction

The remarkable success of deep learning frameworks for numerous inference tasks is well established LeCun et al. (2015). Motivated by the practical implications of the gaps between the empirical observations and theoretical foundations of deep learning, many recent works have explored various approaches to rigorously understand the theory of deep learning models. Multi-layer neural networks have been analyzed extensively in the mean-field regime Mei et al. (2018, 2019); Sirignano & Spiliopoulos (2020). The random features model has also been studied to capture the effects of the regime of parameterization and study phenomenon such as generalization, and "double descent" (see e.g., Mei & Montanari (2022);Lin & Dobriban (2021);Adlam & Pennington (2020)). Among such approaches, the NTKs, first introduced in Jacot et al. (2018), have commonly been used to study the behavior of over-parameterized neural networks Cao et al. (2021); Bietti & Mairal (2019); and are informally defined next.

**Neural Tangent Kernel.** For any given predictor \(f(\mathbf{x};\mathbf{w}):\mathbb{R}^{n\times 1}\times\mathbb{R}^{p}\to \mathbb{R}\) the NTK is the kernel matrix \(\Theta\) defined by the gradient of the predictor output, \(f(\mathbf{x};\mathbf{w})\), with respect to its learnable parameters, \(\mathbf{w}\), as

\[\Theta_{(\mathbf{x}_{i},\mathbf{x}_{j})}(\mathbf{w}):=\langle\nabla_{\mathbf{ w}}f(\mathbf{x}_{i};\mathbf{w}),\nabla_{\mathbf{w}}f(\mathbf{x}_{j};\mathbf{w}) \rangle\, \tag{1}\]

where \(f(\mathbf{x};\mathbf{w})\) represents the predictor output for input data point \(\mathbf{x}\in\mathbb{R}^{n\times 1}\) with the learnable parameters represented by \(\mathbf{w}\in\mathbb{R}^{p}\). The typical setting to study NTKs is that of neural networks in the asymptote of infinite width, where the NTK is constant with respect to the learnable parameters
during training, in contrast to the finite-width scenario Jacot et al. (2018). This constancy of the NTK is a result of certain neural networks transitioning to linearity as their width goes to infinity Liu et al. (2020). NTKs have been leveraged to gain theoretical insights on the behavior of neural networks such as over-parameterized neural networks achieving zero training loss over a non-convex loss landscape Du et al. (2018), the spectral bias of neural networks Cao et al. (2021) and the inductive biases of different neural network architectures Bietti & Mairal (2019).

In particular, the eigenspectrum of the NTK kernel has been linked with the convergence rate of gradient descent for an over-parameterized deep learning model Liu et al. (2022); Arora et al. (2019); Wang et al. (2022a). For instance, gradient descent can achieve faster convergence for a supervised learning problem if the vector of output labels, \(\mathbf{y}\), aligns well with the dominant eigenvectors of the NTK matrix \(\Theta\) Arora et al. (2019); Wang et al. (2022a). For the regression problem pertaining to predicting \(\mathbf{y}\) from \(\mathbf{x}\), our analysis in Section 2 demonstrates that

\[\text{Convergence of gradient descent}\propto\mathbf{y}^{\mathsf{T}}\Theta \mathbf{y} \tag{2}\]

By leveraging the observation above as a motivation, we define \(\mathbf{y}^{\mathsf{T}}\Theta\mathbf{y}\) as Alignment \(\mathcal{A}\).

**GNNs and NTKs.** Given that the NTK \(\Theta\) depends on input data \(\mathbf{x}\) (see equation 1), the alignment \(\mathcal{A}\) inherently captures some version of covariance between output \(\mathbf{y}\) and input \(\mathbf{x}\). Thus, if the NTK \(\Theta\) is'structured' for a given predictor \(f\), the alignment \(\mathcal{A}\) could be leveraged to provide further insights into the design of the predictor \(f\). GNNs are an example of a class of predictors for which the NTK is a function of the graph representation or the graph shift operator (GSO) and input data \(\mathbf{x}\) Krishnagopal & Ruiz (2023). GNNs operating on covariance matrices derived from the input data have been studied previously in Sihag et al. (2022), albeit without any consideration of the insights that could be drawn from the NTKs regarding the choice of graph derived from the data for a supervised learning problem. Many of the existing works that analyze NTKs for GNNs focus on explaining empirically observed trends for GNNs (see Appendix B for expanded literature review).

In this paper, we leverage the structure of NTKs for GNNs to theoretically motivate the choice of a particular GSO. Specifically, if the NTK \(\Theta\) for a GNN is considered to be a function of the form \(\Theta(S,\mathbf{x})\) for a GSO \(S\), the alignment can be represented as \(\mathcal{A}(S,\mathbf{x},\mathbf{y})\), i.e., as a function of the input data \(\mathbf{x}\), output data \(\mathbf{y}\) and \(S\). It is then apparent that optimizing the alignment \(\mathcal{A}\) for a GNN can inform the choice of the GSO \(S\) for a given dataset. A key observation made in this paper is that the alignment \(\mathcal{A}\) is characterized by the cross-covariance between the input and the output and as a result, the optimal GSO for statistical inference is closely related to the cross-covariance.

**Contributions.** In this paper, we consider the setting where the predictor \(f\) is a GNN with graph filter as the convolution operator Gama et al. (2020). Our theoretical contributions in this context are summarized next.

* Our analysis of alignment \(\mathcal{A}\) with graph filter as the predictor motivates cross-covariance between the input and output data as the graph. More precisely, we pose an optimization problem with alignment \(\mathcal{A}\) as the objective function and demonstrate that using the cross-covariance as the GSO maximizes a lower bound on this objective.
* We further extend the results from the graph filter to the scenario of a two-layer GNN as the predictor. Our results show that under certain assumptions, the cross-covariance between the input and the oupt optimizes a lower bound on the alignment for the GNN that has \(\tanh\) activation function. Thus, our analysis motivates using cross-covariance based graphs for a GNN as well.

We validated the insights drawn fron our theoretical results via experiments on the publicly available resting state functional magnetic resonance imaging (rfMRI) data from the Human Connectome Project-Young Adult (HCP-YA) dataset Van Essen et al. (2012). In particular, we considered the task of time series prediction and observed that the GNNs that operated on the cross-covariance between the input and output data achieved better convergence and generalization than those that used the covariance matrix only from the input data.

## 2 Alignment and Convergence of Gradient Descent

In this section, we formalize the concept of alignment \(\mathcal{A}\) and demonstrate its relationship with the convergence of gradient descent for a regression problem. Consider a dataset \(\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{M}\), where \(\mathbf{x}_{i}\in\mathbb{R}^{n\times 1}\), \(\mathbf{y}_{i}\in\mathbb{R}^{n\times 1}\). We aim to leverage the inputs \(\mathbf{x}_{i}\) to estimate the outputs \(\mathbf{y}_{i}\) using a
predictor denoted by \(\mathbf{f}:\mathbb{R}^{n\times 1}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}\). We use the notation \(\mathbf{h}\in\mathbb{R}^{p}\) to denote the vector of all learnable parameters of the predictor. To emphasize the dependence of the predictor \(\mathbf{f}\) on the parameters \(\mathbf{h}\), we use the notation \(\mathbf{f}_{\mathbf{x}_{i}}(\mathbf{h})\) for \(\mathbf{f}(\mathbf{x}_{i},\mathbf{h})\) subsequently. Also the parameters are initialized randomly from a Gaussian distribution \(\mathbf{h}^{(0)}\sim\mathcal{N}(0,\kappa^{2}I)\), where the constant \(\kappa\) controls the magnitude of the initialized parameters. The objective is to minimize the mean squared error (MSE) loss function, defined as

\[\Phi(\mathbf{h})\triangleq\min_{\mathbf{h}\in\mathbb{R}^{p}}\frac{1}{2}\sum_{i=1}^{M} ||\mathbf{y}_{i}-\mathbf{f}_{\mathbf{x}_{i}}(\mathbf{h})||_{2}^{2}. \tag{3}\]

For this purpose, we consider a gradient descent based optimization framework with a learning rate \(\eta>0\). The evolution of the predictor output for a single input \(\mathbf{x}_{i}\) is given by

\[\mathbf{f}_{\mathbf{x}_{i}}(\mathbf{h}^{(t+1)})=\mathbf{f}_{\mathbf{x}_{i}}(\mathbf{h}^{(t)}- \eta\cdot\nabla\Phi(\mathbf{h}^{(t)})) \tag{4}\]

where \(t\) denotes the \(t\)-th step or epoch of gradient descent. To characterize the evolution of the predictor output over the entire dataset, we provide the following definitions

\[\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h})\triangleq\left[\left[\mathbf{f}_{\mathbf{x}_ {1}}(\mathbf{h})\right]^{\mathsf{T}},\,\left[\mathbf{f}_{\mathbf{x}_{2}}(\mathbf{h}) \right]^{\mathsf{T}},\cdots,\,\left[\mathbf{f}_{\mathbf{x}_{M}}(\mathbf{h})\right]^{ \mathsf{T}}\right]^{\mathsf{T}}, \tag{5}\]

\[\tilde{\mathbf{x}}\triangleq\left[\mathbf{x}_{1}^{\mathsf{T}},\,\,\mathbf{x}_ {2}^{\mathsf{T}},\,\cdots,\,\mathbf{x}_{M}^{\mathsf{T}}\right]^{\mathsf{T}},\, \,\,\tilde{\mathbf{y}}\triangleq\left[\mathbf{y}_{1}^{\mathsf{T}},\,\,\mathbf{ y}_{2}^{\mathsf{T}},\,\cdots,\,\mathbf{y}_{M}^{\mathsf{T}}\right]^{\mathsf{T}} \tag{6}\]

where \(\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h}),\tilde{\mathbf{x}}\), and \(\tilde{\mathbf{y}}\) are vectors of length \(nM\). We also define the NTK matrix \(\tilde{\mathbf{\Theta}}(\mathbf{h})\in\mathbb{R}^{nM\times nM}\), which consists of \(M^{2}\) number of \(n\times n\) blocks, such that, the \((i,j)\)-th block is the matrix \(\Theta(\mathbf{x}_{i},\mathbf{x}_{j})\in\mathbb{R}^{n\times n}\) and is given by

\[\Theta(\mathbf{x}_{i},\mathbf{x}_{j})\triangleq\mathrm{J}_{\mathbf{f}_{\mathbf{x} _{i}}}(\mathbf{h}^{(t)})\big{(}\mathrm{J}_{\mathbf{f}_{\mathbf{x}_{j}}}(\mathbf{h}^{(t)}) \big{)}^{\mathsf{T}}. \tag{7}\]

In equation 7, \(\mathrm{J}_{\mathbf{f}_{\mathbf{x}_{i}}}(\mathbf{h})\) denotes the Jacobian matrix with its \((a,b)\)-th entry being \((\mathrm{J}_{\mathbf{f}_{\mathbf{x}_{i}}}(\mathbf{h}))_{ab}=\frac{\partial(\mathbf{f}_{ \mathbf{x}_{i}}(\mathbf{h}))_{a}}{\partial h_{b}}\). If the step size \(\eta\) from equation 4 is sufficiently small, the function \(\mathbf{f}_{\mathbf{x}_{i}}(\mathbf{h}^{(t)})\) can be linearized at each step. In this scenario, the linearized version of the evolution of the predictor output in equation 4 is

\[\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h}^{(t+1)})=\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h} ^{(t)})-\eta\cdot\tilde{\mathbf{\Theta}}(\mathbf{h}^{(t)})\cdot(\tilde{\mathbf{f}}_{ \mathbf{x}}(\mathbf{h}^{(t)})-\tilde{\mathbf{y}}). \tag{8}\]

A typical setting of interest in the existing literature is that of the NTK \(\tilde{\mathbf{\Theta}}(\mathbf{h}^{(t)})\) being a constant with respect to \(\mathbf{h}^{(t)}\). This is because the NTK converges to a constant for many neural networks in the infinite width limit Liu et al. (2020). Theorem 1 characterizes the convergence of gradient descent for the considered multivariate regression problem in this setting (also see Arora et al. (2019), Wang et al. (2022a)). The NTK that is constant with respect to \(\mathbf{h}^{(t)}\) is denoted by \(\tilde{\mathbf{\Theta}}\).

**Theorem 1**.: _In the multivariate regression setting, as described in the beginning of section 2, if the NTK \(\tilde{\mathbf{\Theta}}(\mathbf{h}^{(t)})\) is constant during training and \(\kappa=\mathcal{O}(\varepsilon\sqrt{\frac{\delta}{nM}})\), then with probability at least \(1-\delta\), the training error after \(t\) steps of gradient descent is bounded as_

\[\tilde{\mathbf{y}}^{\mathsf{T}}\left(I-2t\eta\cdot\tilde{\mathbf{\Theta}}\right) \tilde{\mathbf{y}}\pm\mathcal{O}(\varepsilon)\leq||\tilde{\mathbf{f}}_{\mathbf{x} }(\mathbf{h}^{(t)})-\tilde{\mathbf{y}}||_{2}^{2}\leq\tilde{\mathbf{y}}^{\mathsf{T} }\left(I-\eta\cdot\tilde{\mathbf{\Theta}}\right)\tilde{\mathbf{y}}\pm\mathcal{O}(\varepsilon)\]

**Remark 1**.: _In this paper, we primarily consider two classes of predictors. The first class is that of a linear predictor, for which the NTK is a constant given the definition in equation 7. The second class of predictors is that of infinitely wide neural networks (GNNs in particular). We refer the reader to Appendix F and Liu et al. (2020) for a detailed discussion of when and why the NTK can be a constant for neural networks._

Since the term \(\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{\Theta}}\tilde{\mathbf{y}}\) characterizes the upper and lower bounds, the loss \(||\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h}^{(t)})-\tilde{\mathbf{y}}||_{2}^{2}\) is proportional to this term. Based on Theorem 1, we formalize the alignment in Definition 1. A similar definition can be found in Wang et al. (2022a) in the context of active learning.

**Definition 1** (Alignment).: _The alignment between the output \(\tilde{\mathbf{y}}\) and NTK \(\tilde{\mathbf{\Theta}}\) is defined as_

\[\mathcal{A}\triangleq\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{\Theta}}\tilde{ \mathbf{y}}\]
The alignment \(\mathcal{A}\) can be perceived as a metric of correlation between output data and the NTK, and is a characteristic of learning with gradient descent. Using Definition 1, Theorem 1 can be restated as

\[\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{y}}-2t\eta\cdot\mathcal{A}\pm \mathcal{O}(\varepsilon)\leq||\tilde{\mathbf{f}}_{\mathbf{x}}(\mathbf{h}^{(t)})- \tilde{\mathbf{y}}||_{2}^{2}\leq\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{ y}}-\eta\cdot\mathcal{A}\pm\mathcal{O}(\varepsilon) \tag{9}\]

Equation (9) shows that the convergence of gradient descent is positively correlated with \(\mathcal{A}\).

Recall that the NTK \(\tilde{\mathbf{\Theta}}\) is a function of the input data \(\tilde{\mathbf{x}}\) and the learning model \(\mathbf{f}\), even when constant with respect to \(\mathbf{h}^{(t)}\). Therefore, maximizing \(\mathcal{A}\) is contingent on maximizing some kind of cross-covariance between the output data \(\tilde{\mathbf{y}}\) and a function of the input data \(\tilde{\mathbf{x}}\), where the function depends by the learning model \(\mathbf{f}\). This observation motivates us to study the setting where the predictor \(\mathbf{f}\) is a GNN, as a GNN architecture can provide appropriate structure to analyze the connection between alignment, cross-covariance and the structure of the network.

## 3 Alignment motivates Cross-Covariance in GNNs

In this paper, we consider the GNNs for which the convolution operation is a graph filter. A graph filter is characterized by a linear-shift-and-sum operation on the input data and is representative of a large family of convolution operations in GNNs (see the section 'implementation of GCNNs' from Gama et al. (2020)). We begin with the setting where \(\mathbf{f}_{\mathbf{x}}(\mathbf{h})\) is a graph filter.

### NTK and alignment for Graph Filter Model

The formal definition for a graph filter is provided in Definition 1 Gama et al. (2020).

**Definition 2** (Graph Filter).: _Consider a symmetric GSO \(S\in\mathbb{R}^{n\times n}\). A graph filter processes an input \(\mathbf{x}\in\mathbb{R}^{n}\) via a linear-shift-and-sum operation characterized by \(S\), such that, its output is_

\[\mathbf{f}_{\mathbf{x}}(\mathbf{h})=\sum_{k=0}^{K-1}h_{k}S^{k}\mathbf{x}=H(S)\mathbf{ x}\,\quad\text{where}\quad H(S)\triangleq\sum_{k=0}^{K-1}h_{k}S^{k}\, \tag{10}\]

_and \(\mathbf{h}=\{h_{0},h_{1},\cdots,h_{K-1}\}\) is the set of scalars, also referred to as the filter taps or coefficients._

Recall from equation 7 that \(\tilde{\mathbf{\Theta}}(\mathbf{h}^{(t)})\) is a function of the Jacobian matrix \(\mathbf{J}_{\mathbf{f}_{\mathbf{\kappa}_{i}}}(\mathbf{h}^{(t)})\), given by

\[\mathbf{J}_{\mathbf{f}_{\mathbf{\kappa}_{i}}}(\mathbf{h}^{(t)})=\left[\mathbf{x}_{i}|S \mathbf{x}_{i}|S^{2}\mathbf{x}_{i}|\cdots|S^{K-1}\mathbf{x}_{i}\right]. \tag{11}\]

Using equation 11, for any pair of input vectors \((\mathbf{x}_{i},\mathbf{x}_{j})\), the \((i,j)\)-the block of the NTK \(\tilde{\mathbf{\Theta}}(\mathbf{h}^{(t)})\) for a graph filter is given by \(\Theta_{(\mathbf{x}_{i},\mathbf{x}_{j})}(\mathbf{h}^{(t)})=\sum_{k=0}^{K-1}S^{k} \mathbf{x}_{i}(S^{k}\mathbf{x}_{j})^{\mathsf{T}}\). Since the graph filter is a linear model, \(\Theta_{(\mathbf{x}_{i},\mathbf{x}_{j})}(\mathbf{h}^{(t)})\) is independent of \(\mathbf{h}^{(t)}\). Consequently, the NTK \(\Theta_{(\mathbf{x}_{i},\mathbf{x}_{j})}(\mathbf{h}^{(t)})\) for a graph filter is a constant with respect to \(\mathbf{h}^{(t)}\). Next, we provide the NTK for a graph filter.

**Proposition 1** (NTK for a graph filter).: _The NTK for a graph filter is given by_

\[\tilde{\mathbf{\Theta}}_{\text{fit}}(\mathbf{h}^{(t)})=\sum_{k=0}^{K-1}\tilde{S}^{k} \tilde{\mathbf{x}}\tilde{\mathbf{x}}^{\mathsf{T}}\tilde{S}^{k}. \tag{12}\]

_where \(\tilde{S}\in\mathbb{R}^{nM\times nM}\) is a block diagonal matrix consisting of \(M\) blocks of matrix \(S\) on the diagonal and zeros everywhere else._

Given equation 12, we further investigate the impact of shift operator \(S\) on the alignment \(\mathcal{A}\). Also we define the data matrices \(X,Y\) where \(X\) is the input data matrix where the \(i\)-th column is equal to \(\mathbf{x}_{i}\) and similarly for \(Y\). From equation 12, note that the NTK is independent of the filter coefficients \(\mathbf{h}\). As a consequence, \(\mathcal{A}\) for a graph filter (denoted by \(\mathcal{A}_{\text{fit}}\)) depends on the shift operator \(S\) and dataset \((X,Y)\) as follows

\[\mathcal{A}_{\text{fit}}(S,X,Y)=\tilde{\mathbf{y}}^{\mathsf{T}}\left(\sum_{k=0} ^{K-1}\tilde{S}^{k}\tilde{\mathbf{x}}\tilde{\mathbf{x}}^{\mathsf{T}}\tilde{S}^ {k}\right)\tilde{\mathbf{y}}=\sum_{k=0}^{K-1}\left(\tilde{\mathbf{y}}^{\mathsf{ T}}\tilde{S}^{k}\tilde{\mathbf{x}}\right)^{2}=\sum_{k=0}^{K-1}\left(\mathbf{ tr}(Y^{\mathsf{T}}S^{k}X)\right)^{2} \tag{13}\]
The equivalence between different terms in equation 13 follows from the symmetry of \(\tilde{S}\) and the fact that \(\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{S}^{k}\tilde{\mathbf{x}}\) is a scalar. Since a larger \(\mathcal{A}_{\text{fit}}\) is correlated with faster convergence of gradient descent (see equation 9), we further investigate whether the alignment \(\mathcal{A}_{\text{fit}}\) can be optimized by appropriate selection of shift operator matrix \(S\). The objective to optimize \(\mathcal{A}_{\text{fit}}\) can be stated as follows.

\[S^{*}=\operatorname*{arg\,max}_{S}\;\sum_{k=0}^{K-1}\left(\tilde{\mathbf{y}}^{ \mathsf{T}}\tilde{S}^{k}\tilde{\mathbf{x}}\right)^{2}\;\;\text{s.t.}\;\eta\;.\;||\tilde{\Theta}_{\text{fit}}||_{\text{op}}<\alpha\;. \tag{14}\]

The constraint \(||\eta\cdot\tilde{\Theta}_{\text{fit}}||_{\text{op}}<\alpha\), for some \(\alpha>0\), in equation 14 is necessary to ensure the convergence of gradient descent. This constraint also eliminates trivial solutions (such as multiplying a given \(S\) with an arbitrarily large positive constant to inflate \(\mathcal{A}\) in isolation). The optimization problem in equation 14, while meaningful, can be analytically intractable due to complications arising from the polynomial functions of \(\tilde{S}\) and the objective function and the constraint being non-convex. In order to provide an analytically tractable solution to \(S\), we consider the lower bound on \(\mathcal{A}\) next.

**Lemma 1**.: _[Lower bound on \(\mathcal{A}_{\text{fit}}\).] The alignment \(\mathcal{A}_{\text{fit}}\) satisfies \(\mathcal{A}_{\text{fit}}(S,X,Y)\geq\mathcal{A}_{\text{L}}(S,X,Y)\), where_

\[\mathcal{A}_{\text{L}}(S,X,Y)\triangleq\left(\frac{1}{\sqrt{K}}\boldsymbol{tr }\Big{(}\Big{(}\sum_{k=0}^{K-1}S^{k}\Big{)}C_{XY}\Big{)}\right)^{2},\quad\text {and}\quad C_{XY}\triangleq\frac{1}{2}(XY^{\mathsf{T}}+YX^{\mathsf{T}})\;. \tag{15}\]

Henceforth, we focus on characterizing \(S\) that maximizes \(\mathcal{A}_{\text{L}}(S,X,Y)\). Our experiments in Section 4 also demonstrate that the insights drawn from optimizing \(\mathcal{A}_{\text{L}}\) are practically meaningful. Next, we provide a constraint that depends on the choice of GSO and not on the input data.

**Lemma 2**.: _If the degree \(K\) polynomial in the shift operator \(S\) has a bounded Frobenius norm, the operator norm of the NTK matrix is also bounded as follows:_

\[\left\|\sum_{k=0}^{K-1}S^{k}\right\|_{F}\leq\sqrt{\alpha/(\eta M)}\;\;\Rightarrow \;\;\eta\cdot\left\|\sum_{k=0}^{K-1}\tilde{S}^{k}\tilde{\mathbf{x}}\tilde{ \mathbf{x}}^{\mathsf{T}}\tilde{S}^{k}\right\|_{\text{op}}\leq\alpha \tag{16}\]

The constraint on the left in equation 16 is more straightforward to work with in the analysis since it only depends on \(S\), while also ensuring that the constraint in equation 14 is satisfied. Putting together \(\mathcal{A}_{\text{L}}(S,X,Y)\) and the revised constraint, we get the following _optimization problem_.

\[S^{*}=\operatorname*{arg\,max}_{S}\;\mathcal{A}_{\text{L}}(S,X,Y)\;\;\text{s. t.}\;\left\|\sum_{k=0}^{K-1}S^{k}\right\|_{F}\leq\sqrt{\alpha/(\eta M)}\;. \tag{17}\]

**Theorem 2** (GSO in graph filter.).: _A GSO \(S^{*}\) that satisfies_

\[\sum_{k=0}^{K-1}(S^{*})^{k}=\mu\cdot C_{XY}\;,\quad\text{where}\quad\mu=\frac {\sqrt{\alpha/(\eta M)}}{||C_{XY}||_{F}}\;. \tag{18}\]

_is the solution to the optimization problem in equation 17._

Theorem 2 clearly demonstrates the association between the optimal GSO that optimizes \(\mathcal{A}_{\text{L}}(S,X,Y)\) and \(C_{XY}\), which is a measure of cross-covariance. For instance, if \(K=2\), then it can be concluded from equation 18 that

\[I+S^{*}=\mu\cdot C_{XY}=S^{*}=\mu\cdot C_{XY}-I \tag{19}\]

The observation in equation 19 motivates the potential choice of a normalized cross-covariance matrix as a GSO when graph filter is deployed as the predictor \(\boldsymbol{f}_{\mathbf{x}}(\boldsymbol{h})\). Next, we discuss how this observation extends to the setting where \(\boldsymbol{f}_{\mathbf{x}}(\boldsymbol{h})\) is a GNN.

### NTK and Alignment for GNN

To start with, we formalize the GNN architecture that is the focus of our analysis. The ability to learn non-linear mappings by GNNs is fundamentally based on concatenating an element-wise non-linearity with a graph filter to form a graph perceptron, which is realized via a point-wise
non-linearity \(\sigma(\cdot)\) as \(\sigma(H(S)\mathbf{x})\). In the remainder of this paper, we will focus on a two-layer GNN that admits a single input feature \(\mathbf{x}\in\mathbb{R}^{n}\) and the GNN output is a vector of length \(n\), as dictated by the problem definition in Section 2. The general definition for a GNN, along with additional experimental results for GNNs with depth larger than two layers, have been provided in Appendix H.

**Two-Layer GNN Architecture** (see Fig. 3). In the first layer, the input vector \(\mathbf{x}\in\mathbb{R}^{n}\), is processed by \(F\) graph perceptrons to output \(F\)\(n\)-dimensional outputs given by \(\mathbf{q}_{(1)}^{f},\forall f\in\{1,\cdots,F\}\), as follows

\[\mathbf{u}_{(1)}^{f}=H_{(1)}^{f}(S)\mathbf{x}=\sum_{k=0}^{K-1}h_{(1),k}^{f}S^{ k}\mathbf{x},\forall f\in\{1,\ldots,F\};\quad\text{and}\quad\mathbf{q}_{(1)}^{f}= \sigma\Big{(}\mathbf{u}_{(1)}^{f}\Big{)}. \tag{20}\]

In the second layer, each of the outputs of the previous layer, \(\mathbf{q}_{(1)}^{f}\) are processed by a graph filter as

\[\mathbf{u}_{(2)}^{f}=H_{(2)}^{f}(S)\mathbf{q}_{(1)}^{f}=\sum_{k=0}^{K-1}h_{(2),k}^{f}S^{k}\mathbf{q}_{(1)}^{f},\forall f=\in\{1,\ldots,F\}. \tag{21}\]

Finally, the terms \(\mathbf{u}_{(2)}^{f}\) are aggregated to get the output at the second layer (also the GNN output) as

\[\mathbf{f}_{\mathbf{x}}(\mathbf{h})=\frac{1}{\sqrt{F}}\sum_{f=1}^{F}\mathbf{u}_{(2)} ^{f} \tag{22}\]

The absence of a non-linearity in the final layer (equation 22) is necessary for having a constant NTK in the infinite width limit (see Liu et al. (2020)).

**Proposition 2** (NTK for a two-layer GNN).: _The NTK for the two-layer GNN is given by_

\[\tilde{\mathbf{\Theta}}_{GNN}(\mathbf{h})=\frac{1}{F}\sum_{f=1}^{F}\sum_{k=0}^{K-1} \left(\mathbf{c}_{f,k}^{(1)}\right)\left(\mathbf{c}_{f,k}^{(1)}\right)^{ \mathsf{T}}+\frac{1}{F}\sum_{f=1}^{F}\sum_{k=0}^{K-1}\left(\mathbf{c}_{f,k}^{ (2)}\right)\left(\mathbf{c}_{f,k}^{(2)}\right)^{\mathsf{T}} \tag{23}\]

\[\text{where}\quad\mathbf{c}_{f,k}^{(1)}\triangleq H_{f}(\tilde{S})\cdot\text{diag }(\sigma^{\prime}(G_{f}(\tilde{S})\tilde{\mathbf{x}})\tilde{S}^{k}\tilde{ \mathbf{x}}),\quad\text{and}\quad\mathbf{c}_{f,k}^{(2)}\triangleq\tilde{S}^{k }\sigma\left(G_{f}(\tilde{S})\tilde{\mathbf{x}}\right). \tag{24}\]

_In equation 24, \(\mathbf{c}_{f,k}^{(\ell)}\in\mathbb{R}^{M\times 1}\) is the vector determined by picking out the column that pertains to the derivative of the network output with regards to the parameter indexed by \((f,k,\ell)\), namely, the \(k\)-th coefficient of the \(f\)-th filter in layer \(\ell\), from every Jacobian matrices \(J_{\mathbf{f}_{\mathbf{x}_{i}}},\forall i\in\{1,\cdots M\}\) and stacking all these vectors together._

The NTK in equation 23 is an aggregation of two terms, where the first term is associated with the first layer and the second term with the second layer. It follows from Definition 1 and equation 23 that the alignment for a two-layer GNN is also composed of two terms that represent the two layers. Henceforth, we focus on the results pertaining to the second term in equation 23. This implies that our results correspond to a two-layer GNN where only the parameters of the second layer are trained and the parameters of the first layer are fixed. The analysis (and results) if we also consider the first layer in this analysis is similar and has been relegated to Appendix G.

In the subsequent discussions, the notation \(\tilde{\mathbf{\Theta}}_{GNN}\) denotes the second term in equation 23 when the width of hidden layer approaches infinity i.e., \(F\rightarrow\infty\). Therefore, \(\tilde{\mathbf{\Theta}}_{GNN}\) is given by

\[\tilde{\mathbf{\Theta}}_{GNN}(\mathbf{h}) =\lim_{F\rightarrow\infty}\frac{1}{F}\sum_{f=1}^{F}\sum_{k=0}^{K-1 }\left(\tilde{S}^{k}\sigma\left(G_{f}(\tilde{S})\tilde{\mathbf{x}}\right) \right)\left(\tilde{S}^{k}\sigma\left(G_{f}(\tilde{S})\tilde{\mathbf{x}} \right)\right)^{\mathsf{T}}\] \[=\sum_{k=0}^{K-1}\tilde{S}^{k}\mathop{\mathbb{E}}_{\mathbf{g}\sim \mathcal{N}(0,I)}\left[\sigma\left(G(\tilde{S})\tilde{\mathbf{x}}\right)\left( \sigma\left(G(\tilde{S})\tilde{\mathbf{x}}\right)\right)^{\mathsf{T}}\right] \tilde{S}^{k}=\sum_{k=0}^{K-1}\tilde{S}^{k}E\tilde{S}^{k} \tag{25}\]

The expectation matrix \(E\triangleq\mathop{\mathbb{E}}_{\mathbf{g}\sim\mathcal{N}(0,I)}\left[\sigma\left(G (\tilde{S})\tilde{\mathbf{x}}\right)\left(\sigma\left(G(\tilde{S})\tilde{ \mathbf{x}}\right)\right)^{\mathsf{T}}\right]\) is instrumental for the analysis of the alignment. Before proceeding, we provide the following remark pertinent to the analysis.
**Remark 2**.: _As a byproduct of the output layer being linear, the NTK \(\tilde{\mathbf{\Theta}}_{GNN}\) in equation 25 does not depend on the parameters of the second layer, i.e., \(\mathbf{h}_{f},\forall f\in\{1,\cdots,F\}\). Hence, the NTK in equation 25 could be considered a constant if only the second layer of GNN is trained. For completeness, our discussion in Appendix F demonstrates further that as \(F\to\infty\), the NTK in equation 23 also approaches a constant behavior._

From equation 25, the alignment can be written in terms of \(E\) as follows

\[\mathcal{A}=\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{\Theta}}_{GNN} \tilde{\mathbf{y}}=\tilde{\mathbf{y}}^{\mathsf{T}}\Big{(}\sum_{k=0}^{K-1} \tilde{S}^{k}E\tilde{S}^{k}\Big{)}\tilde{\mathbf{y}}=\mathbf{tr}(QE)\, \tag{26}\]

where we have defined the matrix \(Q\) as \(Q\triangleq\sum_{k=0}^{K-1}\tilde{S}^{k}\tilde{\mathbf{y}}\tilde{\mathbf{y}}^ {\mathsf{T}}\tilde{S}^{k}\). Above, we used the cyclic property of the trace and the fact that \(\tilde{\mathbf{y}}^{\mathsf{T}}\tilde{\mathbf{\Theta}}_{GNN}\tilde{\mathbf{y}}\) is a scalar. In order to evaluate \(E\), we define the vectors \(\mathbf{z}^{(\ell)}\in\mathbb{R}^{K\times 1}\) as \(\mathbf{z}^{(\ell)}\triangleq\Big{[}\tilde{\mathbf{x}}_{\ell},\ (\tilde{S}\tilde{ \mathbf{x}})_{\ell},\ \cdots,\ (\tilde{S}^{K-1}\tilde{\mathbf{x}})_{\ell}\Big{]}^{ \mathsf{T}}\) where \(\ell\in\{1,\cdots,nM\}\) and \((\tilde{S}^{k}\tilde{\mathbf{x}})_{\ell}\) denotes the \(\ell\)-th entry of the vector \(\tilde{S}^{k}\tilde{\mathbf{x}}\). Thus, the \((a,b)\)-th entry of \(E\), i.e., \(E_{ab}\), from equation 25 is

\[E_{ab}=\underset{\mathbf{g}\sim\mathcal{N}(0,I)}{\mathbb{E}}\Big{[}\sigma\Big{(} \langle\mathbf{g},\mathbf{z}^{(a)}\rangle\Big{)}\cdot\sigma\Big{(}\langle\mathbf{g},\mathbf{z }^{(b)}\rangle\Big{)}\Big{]}. \tag{27}\]

**Linear GNNs.** We next discuss the scenario when the function \(\sigma(\cdot)\) is an identity function, i.e., \(\sigma(z)=z\). The results drawn from this setting will be leveraged later in the setting when \(\sigma(\cdot)\) is not an identity function. When \(\sigma(z)=z\), equation 27 reduces to

\[E_{ab}=\underset{\mathbf{g}\sim\mathcal{N}(0,I)}{\mathbb{E}}\Big{[}\langle\mathbf{g}, \mathbf{z}^{(a)}\rangle\cdot\langle\mathbf{g},\mathbf{z}^{(b)}\rangle\Big{]}=\langle\mathbf{z }^{(a)},\mathbf{z}^{(b)}\rangle \tag{28}\]

We denote the matrix \(E\) in this linear setting by \(B_{\text{lin}}\), which is given by

\[(B_{\text{lin}})_{ab}\triangleq\langle\mathbf{z}^{(a)},\mathbf{z}^{(b)}\rangle\Rightarrow B _{\text{lin}}=\sum_{k=0}^{K-1}\tilde{S}^{k}\tilde{\mathbf{x}}\tilde{\mathbf{ x}}^{\mathsf{T}}\tilde{S}^{k} \tag{29}\]

Thus, the alignment in the linear setting is given by

\[\mathcal{A}_{\text{lin}}\triangleq\mathbf{tr}(QB_{\text{lin}})=\sum_{k=0}^{K -1}\tilde{\mathbf{y}}^{T}\tilde{S}^{k}B\tilde{S}^{k}\tilde{\mathbf{y}}=\sum_{k =0}^{K-1}\sum_{k^{\prime}=0}^{K-1}\tilde{\mathbf{y}}^{T}\tilde{S}^{k+k^{ \prime}}\tilde{\mathbf{x}}\tilde{T}^{k}\tilde{S}^{k+k^{\prime}}\tilde{ \mathbf{y}} \tag{30}\]

The analysis of alignment \(\mathcal{A}_{\text{lin}}\) in equation 30 using similar arguments as that for a graph filter in Section 3.1 yields a similar condition on the GSO \(S\) as in Theorem 2. The corollaries provided next formalize this observation. First, the following corollary provides a lower bound on \(\mathcal{A}_{\text{lin}}\).

**Corollary 1**.: _[Lower bound on \(\mathcal{A}_{\text{lin}}\)] The term \(\mathcal{A}_{\text{lin}}=\textbf{tr}(QB_{\text{lin}})\) satisfies \(\mathcal{A}_{\text{lin}}\geq\mathcal{A}_{\mathsf{U}^{\prime}}(S,X,Y)\),_

\[\text{where}\quad\mathcal{A}_{\mathsf{U}^{\prime}}(S,X,Y)\triangleq\Big{(} \frac{1}{\sqrt{K}}\textbf{tr}\Big{(}\Big{(}\sum_{k=0}^{K-1}\sum_{k^{\prime}= 0}^{K-1}(S^{*})^{k+k^{\prime}}\Big{)}C_{XY}\Big{)}\Big{)}^{2}\, \tag{31}\]

Next, we present an _optimization problem_ similar to the one for the graph filter in equation 17 next.

\[S^{*}=\underset{S}{\arg\max}\ \mathcal{A}_{\mathsf{U}^{\prime}}(S,X,Y)\ \textbf{ s.t. }\ \ \left\|\sum_{k=0}^{K-1}\sum_{k^{\prime}=0}^{K-1}S^{k+k^{\prime}}\right\|_{F}\leq \sqrt{\alpha/(\eta M)}. \tag{32}\]

The solution to the optimization problem in equation 32 is presented next.

**Corollary 2**.: _[Extension of Theorem 1 to linear GNN] The GSO \(S^{*}\) that solves the optimization problem in equation 32 must satisfy_

\[\sum_{k=0}^{K-1}\sum_{k=0}^{K-1}(S^{*})^{k+k^{\prime}}=\mu\cdot C_{XY}\,\quad \text{where}\quad\mu=\frac{\sqrt{\alpha/(\eta M)}}{||C_{XY}||_{F}}. \tag{33}\]
Corollary 2 establishes that the cross-covariance \(C_{XY}\) is instrumental to optimizing \(\mathcal{A}_{\mathsf{L}^{\prime}}\) for the considered two-layer GNN architecture when \(\sigma(\cdot)\) is an identity function. In general, this observation holds for linear GNNs of any arbitrary depth.

**GNNs with non-linear activation function.** Next, we investigate the conditions under which the observation in Corollary 2 extends to a more general setting, in which \(\sigma(\cdot)\) is not an identity function. We will focus our theoretical analysis on the case where \(\sigma(z)=\text{tanh}(z)\) and from here on \(\mathcal{A}\) will denote the alignment for this case. The experimental results (see Appendix H) show that similar results hold for some other activation functions like ReLU in practice. First, we evaluate the expectation in equation 25. By leveraging the theory of Hermite polynomials 1, the Hermite expansions of \(\sigma\left(\langle\mathbf{g},\mathbf{z}^{(a)}\rangle\right)\) and \(\sigma\left(\langle\mathbf{g},\mathbf{z}^{(b)}\rangle\right)\) enables the expansion of \(E\) and subsequently \(\mathcal{A}\). These expansions are formalized next.

Footnote 1: See the proof of Lemma 3 for an overview of the Hermite polynomials and how we utilized the Hermite expansion.

**Lemma 3** (Expansion of \(E\) and \(\mathcal{A}\)).: _The Hermite expansion of \(E\) is given by \(E=B+\Delta B\), where \(B\in\mathbb{R}^{nM\times nM}\) represents the first non-zero term in the expansion and \(\Delta B\in\mathbb{R}^{nM\times nM}\) includes all the subsequent terms. For the \((a,b)\)-th element of \(B\) and \(\Delta B\), we have_

\[B_{ab}=\alpha_{1}\beta_{1}\cdot\frac{\langle\mathbf{z}^{(a)},\mathbf{z}^{(b)}\rangle} {||\mathbf{z}^{(a)}||_{2}\cdot||\mathbf{z}^{(b)}||_{2}},\ \text{ and }\ (\Delta B)_{ab}=\sum_{i=3,5,\cdots}^{\infty}\alpha_{i}\beta_{i}\cdot\left( \frac{\langle\mathbf{z}^{(a)},\mathbf{z}^{(b)}\rangle}{||\mathbf{z}^{(a)}||_{2}\cdot||\bm {z}^{(b)}||_{2}}\right)^{i}. \tag{34}\]

_Hence, the alignment \(\mathcal{A}\) in equation 26 admits the expansion_

\[\mathcal{A}=\mathbf{tr}(QE)=\mathbf{tr}(QB)+\mathbf{tr}(Q\Delta B). \tag{35}\]

The scalar coefficients \(\alpha_{i},\beta_{i}\) in equation 34 depend on \(||\mathbf{z}^{(a)}||_{2}\) and \(||\mathbf{z}^{(b)}||_{2}\), respectively and the choice of \(\sigma(\cdot)\). The disintegration of the alignment into two terms in equation 35 is useful because the second term can be shown to be relatively small and, using the observation that in the linear case \(E_{ab}=\langle\mathbf{z}^{(a)},\mathbf{z}^{(b)}\rangle\), \(\mathbf{tr}(QB)\) can be related to \(\mathcal{A}_{\text{fin}}\). We next provide two lemmas relevant to this.

**Lemma 4**.: _Given a family of matrices \(S\in\mathbb{S}^{n\times n}\) that have a bounded norm, \(||S||_{\text{op}}\leq\nu\), we have_

\[\mathbf{tr}(QB)\geq\rho\mathcal{A}_{\text{fin}} \tag{36}\]

_where \(\rho\) is a constant that depends on the choice of non-linearity function \(\sigma(\cdot)\)._

The next lemma shows that the elements of the matrix \(\Delta B\) are smaller compared to corresponding elements in \(B\), which implies that the second term in equation 35 can't decrease \(\mathcal{A}\) too much.

**Lemma 5**.: _Each element of \(\Delta B\) has the same sign as the corresponding element in \(B\). Also, the following element-wise inequality holds between the two matrices:_

\[|\Delta B|\leq\beta\cdot|B| \tag{37}\]

_where \(\beta\) is a constant that depends on our choice of non-linearity and is determined from the proof._

Putting the above two lemmas together, we reach the following conclusion about the alignment of the two-layer GNN with \(\sigma(\cdot)\) as \(\tanh\).

**Theorem 3**.: _Given a family of matrices \(S\in\mathbb{S}^{n\times n}\) that have a bounded norm, \(||S||_{\text{op}}\leq\nu\) and that satisfy \(\mathcal{A}_{\text{fin}}=\mathbf{tr}\left(QB_{\text{fin}}\right)\geq\xi\cdot||Q|| _{F}||B_{\text{fin}}||_{F}\) for some constant \(0<\xi\leq 1\), \(\mathcal{A}_{\text{fin}}\) lower bounds the alignment for the two-layer GNN with \(\tanh\) non-linearity, \(\mathcal{A}\), up to a constant as follows_

\[\mathcal{A}\geq\Big{(}c-\frac{d}{\xi}\Big{)}\mathcal{A}_{\text{fin}}\, \tag{38}\]

_for some positive constants \(c\) and \(d\)._

**Remark 3** (GSO in GNNs).: _Our key takeaway from Theorem 3 is that when maximizing \(\mathcal{A}_{\text{fin}}\) over the family of shift operators that satisfy the assumption with a sufficiently large \(\xi\), such that \(c-d/\xi\) is a positive constant, we are essentially maximizing a lower bound on the alignment \(\mathcal{A}\) of the two-layer GNN. This observation, together with Corollary 3, motivates using the cross-covariance matrix \(C_{XY}\) as a GSO for the two-layer GNN._
**Alignment, The NTK and Generalization** Thus far, we have provided the theoretical results motivated by the fact that larger alignment can imply faster convergence of gradient descent during training. However, the alignment and NTK are also closely related to generalization. Specifically, the analyses pertaining to generalization from Arora et al. (2019) and Wang et al. (2022) can be extended to the case of graph filters, which leads to the conclusion that larger alignment could also lead to smaller generalization error. Hence, the results on improved training and generalization together motivate models with larger alignment in practice. The analysis regarding generalization has been provided in Appendix E.

## 4 Experiments

In this section, we provide the experiments that validate the theoretical insights pertaining to the cross-covariance matrix being an optimal GSO for GNN training and generalization with respect to GSO derived only from the input data for a regression task. The dataset and inference task for this purpose are described below.

**Data.** The HCP-YA dataset is a publicly available brain imaging dataset collected over a population of 1003 healthy adults in the age range of 22-35 years Van Essen et al. (2012, 2013). In our experiments, we leveraged the rfMRI data for each subject made available by HCP. This data consisted of a multi-variate time series of \(100\) features, with each time series consisting of 4500 time points.

**Inference task.** Noting that the \(100\) features could be considered as \(100\) nodes of a graph, our objective was to use the data at all nodes at the current time step for an individual to predict the data at all nodes at a future time step. Specifically, given the signal value at time step \(t\) as

Figure 1: Experimental Results
for an individual, we aimed to predict the signal value after \(\Delta t\) time steps, i.e., \(\mathbf{z}^{(t+\Delta t)}\in\mathbb{R}^{100}\). For every \(\Delta t\in\{1,2,3,4,5\}\), a separate training/test set of size \(N_{train}=1000,\ N_{test}=100\) was created, such that, for the signal at a time point \(t\), i.e., \(\mathbf{z}^{(t)}\) as the input, the signal after \(\Delta t\) time steps \(\mathbf{z}^{(t+\Delta t)}\) was the output to be predicted. For additional implementation details, see Appendix H.

**Performance evaluation.** We trained two sets of GNNs and two sets of graph filters using the time series data of each individual for a given \(\Delta t\), where one set comprised of predictors with \(C_{XY}\) as the GSO and the other with \(C_{XX}\) as the GSO. The GNNs with \(C_{XX}\) as the GSO have been studied before as VNNs in Sihag et al. (2022) and provide an appropriate baseline for comparison as it is representative of GNNs with GSOs extracted only from the input data. Figures 0(a) and 0(b) illustrate faster convergence of both training loss and test loss during gradient descent for predictors with \(C_{XY}\) as compared to those with \(C_{XX}\) for one representative individual when \(\Delta t=1\). This observation was consistent for graph filters and GNNs. For each individual, the training process for every architecture was repeated 10 times. The average of these runs is shown in Figures 0(a) and 0(b).

Further, we checked whether these observations were consistent across the dataset and for different \(\Delta t\). Figure 0(c) illustrates the gap between the test error for predictors with \(C_{XY}\) and \(C_{XX}\) and different values of \(\Delta t\), averaged across all individuals. Even as the accuracy of prediction diminished with increasing \(\Delta t\), we observed a consistent gain in test performance when using \(C_{XY}\) as compared to \(C_{XX}\). Similarly, Fig. 0(d) shows that predictors with \(C_{XY}\) achieved smaller training error relative to those with \(C_{XX}\) at each epoch of gradient descent, averaged across the dataset. Thus, the observations in Fig. 1 validate the theoretical insights drawn from the analysis that argued for \(C_{XY}\) as an appropriate GSO for GNNs that can achieve smaller training error and better generalization. We refer the reader to Appendix C for the conclusions, limitations, and potential future directions of our work.

## 5 Reproducibility Statement

**Theoretical Proofs.** The proofs for all theorems and lemmas in the main body of the paper can be found in Appendix D.

**Experiments.** Additional experimental details and results have been provided in Appendix H. In addition, the code and data for training the relevant models and producing results similar to what is shown in Figures 0(a) and 0(b) for one individual from the HCP-YA dataset can be found in [https://github.com/shervinkh2000/Cross_Covariance_NTK](https://github.com/shervinkh2000/Cross_Covariance_NTK). The complete HCP-YA dataset is publicly available and accessible from [https://db.humanconnectome.org/](https://db.humanconnectome.org/) as per the data use terms.

**Additional Details.** Comments made regarding the relation between the NTK and generalization made in the main body of the paper have been thoroughly explained in Appendix E. Additional necessary theoretical considerations that did not directly contribute to the main message of the paper, namely, that the NTK of the GNN architecture that we analyzed is constant in the infinite-width limit, and that training the first layer of the GNN leads to similar results to those we saw for the second layer in section 3.2, have been discussed further in Appendices F and G respectively.
