Graph Neural Network-Inspired Kernels for Gaussian Processes in
  Semi-Supervised Learning
Gaussian processes (GPs) are an attractive class of machine learning models
because of their simplicity and flexibility as building blocks of more complex
Bayesian models. Meanwhile, graph neural networks (GNNs) emerged recently as a
promising class of models for graph-structured data in semi-supervised learning
and beyond. Their competitive performance is often attributed to a proper
capturing of the graph inductive bias. In this work, we introduce this
inductive bias into GPs to improve their predictive performance for
graph-structured data. We show that a prominent example of GNNs, the graph
convolutional network, is equivalent to some GP when its layers are infinitely
wide; and we analyze the kernel universality and the limiting behavior in
depth. We further present a programmable procedure to compose covariance
kernels inspired by this equivalence and derive example kernels corresponding
to several interesting members of the GNN family. We also propose a
computationally efficient approximation of the covariance matrix for scalable
posterior inference with large-scale data. We demonstrate that these
graph-based kernels lead to competitive classification and regression
performance, as well as advantages in computation time, compared with the
respective GNNs.
Zehao Niu, Mihai Anitescu, Jie Chen
# Graph Neural Network-Inspired Kernels for Gaussian Processes in Semi-Supervised Learning

###### Abstract

Gaussian processes (GPs) are an attractive class of machine learning models because of their simplicity and flexibility as building blocks of more complex Bayesian models. Meanwhile, graph neural networks (GNNs) emerged recently as a promising class of models for graph-structured data in semi-supervised learning and beyond. Their competitive performance is often attributed to a proper capturing of the graph inductive bias. In this work, we introduce this inductive bias into GPs to improve their predictive performance for graph-structured data. We show that a prominent example of GNNs, the graph convolutional network, is equivalent to some GP when its layers are infinitely wide; and we analyze the kernel universality and the limiting behavior in depth. We further present a programmable procedure to compose covariance kernels inspired by this equivalence and derive example kernels corresponding to several interesting members of the GNN family. We also propose a computationally efficient approximation of the covariance matrix for scalable posterior inference with large-scale data. We demonstrate that these graph-based kernels lead to competitive classification and regression performance, as well as advantages in computation time, compared with the respective GNNs.

## 1 Introduction

Gaussian processes (GPs) (Rasmussen and Williams, 2006) are widely used in machine learning, uncertainty quantification, and global optimization. In the Bayesian setting, a GP serves as a prior probability distribution over functions, characterized by a mean (often treated as zero for simplicity) and a covariance. Conditioned on observed data with a Gaussian likelihood, the random function admits a posterior distribution that is also Gaussian, whose mean is used for prediction and the variance serves as an uncertainty measure. The closed-form posterior allows for exact Bayesian inference, resulting in great attractiveness and wide usage of GPs.

The success of GPs in practice depends on two factors: the observations (training data) and the covariance kernel. We are interested in semi-supervised learning, where only a small amount of data is labeled while a large amount of unlabeled data can be used together for training (Zhu, 2008). In recent years, graph neural networks (GNNs) (Zhou et al., 2020; Wu et al., 2021) emerged as a promising class of models for this problem, when the labeled and unlabeled data are connected by a graph. The graph structure becomes an important inductive bias that leads to the success of GNNs. This inductive bias inspires us to design a GP model under limited observations, by building the graph structure into the covariance kernel.

An intimate relationship between neural networks and GPs is known: a neural network with fully connected layers, equipped with a prior probability distribution on the weights and biases, converges to a GP when each of its layers is infinitely wide (Lee et al., 2018; de G. Matthews et al., 2018). Such a result is owing to the central limit theorem (Neal, 1994; Williams, 1996) and the GP covariance
can be recursively computed if the weights (and biases) in each layer are iid Gaussian. Similar results for other architectures, such as convolution layers and residual connections, were subsequently established in the literature (Novak et al., 2019; Garriga-Alonso et al., 2019).

One focus of this work is to establish a similar relationship between GNNs and the limiting GPs. We will derive the covariance kernel that incorporates the graph inductive bias as GNNs do. We start with one of the most widely studied GNNs, the graph convolutional network (GCN) (Kipf and Welling, 2017), and analyze the kernel universality as well as the limiting behavior when the depth also tends to infinity. We then derive covariance kernels from other GNNs by using a programmable procedure that corresponds every building block of a neural network to a kernel operation.

Meanwhile, we design efficient computational procedures for posterior inference (i.e., regression and classification). GPs are notoriously difficult to scale because of the cubic complexity with respect to the number of training data. Benchmark graph datasets used by the GNN literature may contain thousands or even millions of labeled nodes (Hu et al., 2020). The semi-supervised setting worsens the scenario, as the covariance matrix needs to be (recursively) evaluated in full because of the graph convolution operation. We propose a Nystrom-like scheme to perform low-rank approximations and apply the approximation recursively on each layer, to yield a low-rank kernel matrix. Such a matrix can be computed scalably. We demonstrate through numerical experiments that the GP posterior inference is much faster than training a GNN and subsequently performing predictions on the test set.

We summarize the contributions of this work as follows:

1. We derive the GP as a limit of the GCN when the layer widths tend to infinity and study the kernel universality and the limiting behavior in depth.
2. We propose a computational procedure to compute a low-rank approximation of the covariance matrix for practical and scalable posterior inference.
3. We present a programmable procedure to compose covariance kernels and their approximations and show examples corresponding to several interesting members of the GNN family.
4. We conduct comprehensive experiments to demonstrate that the GP model performs favorably compared to GNNs in prediction accuracy while being significantly faster in computation.

## 2 Related Work

It has long been observed that GPs are limits of standard neural networks with one hidden layer when the layer width tends to infinity (Neal, 1994; Williams, 1996). Recently, renewed interests in the equivalence between GPs and neural networks were extended to deep neural networks (Lee et al., 2018; de G. Matthews et al., 2018) as well as modern neural network architectures, such as convolution layers (Novak et al., 2019), recurrent networks (Yang, 2019), and residual connections (Garriga-Alonso et al., 2019). The term NNGP (neural network Gaussian process) henceforth emerged under the context of Bayesian deep learning. Besides the fact that an infinite neural network defines a kernel, the training of a neural network by using gradient descent also defines a kernel--the neural tangent kernel (NTK)--that describes the evolution of the network (Jacot et al., 2018; Lee et al., 2019). Library supports in Python were developed to automatically construct the NNGP and NTK kernels based on programming the corresponding neural networks (Novak et al., 2020).

GNNs are neural networks that handle graph-structured data (Zhou et al., 2020; Wu et al., 2021). They are a promising class of models for semi-supervised learning. Many GNNs use the message-passing scheme (Gilmer et al., 2017), where neighborhood information is aggregated to update the representation of the center node. Representative examples include GCN (Kipf and Welling, 2017), GraphSAGE (Hamilton et al., 2017), GAT (Velickovic et al., 2018), and GIN (Xu et al., 2019). It is found that the performance of GNNs degrades as they become deep; one approach to mitigating the problem is to insert residual/skip connections, as done by JumpingKnowledge (Xu et al., 2018), APPNP (Gasteiger et al., 2019), and GCNII (Chen et al., 2020).

GP inference is too costly, because it requires the inverse of the \(N\times N\) dense kernel matrix. Scalable approaches include low-rank methods, such as Nystrom approximation (Drineas and Mahoney, 2005), random features (Rahimi and Recht, 2007), and KISS-GP (Wilson and Nickisch, 2015); as well as multi-resolution (Katzfuss, 2017) and hierarchical methods (Chen et al., 2017; Chen and Stein, 2021).
Prior efforts on integrating graphs into GPs exist. Ng et al. (2018) define a GP kernel by combing a base kernel with the adjacency matrix; it is related to a special case of our kernels where the network has only one layer and the output admits a robust-max likelihood for classification. Hu et al. (2020) explore a similar route to us, by taking the limit of a GCN, but its exploration is less comprehensive because it does not generalize to other GNNs and does not tackle the scalability challenge.

## 3 Graph Convolutional Network as a Gaussian Process

We start with a few notations used throughout this paper. Let an undirected graph be denoted by \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) with \(N=|\mathcal{V}|\) nodes and \(M=|\mathcal{E}|\) edges. For notational simplicity, we use \(A\in\mathbb{R}^{N\times N}\) to denote the original graph adjacency matrix or any modified/normalized version of it. Using \(d_{l}\) to denote the width of the \(l\)-th layer, the layer architecture of GCN reads

\[X^{(l)}=\phi\left(AX^{(l-1)}W^{(l)}+b^{(l)}\right), \tag{1}\]

where \(X^{(l-1)}\in\mathbb{R}^{N\times d_{l-1}}\) and \(X^{(l)}\in\mathbb{R}^{N\times d_{l}}\) are layer inputs and outputs, respectively; \(W^{(l)}\in\mathbb{R}^{d_{l-1}\times d_{l}}\) and \(b^{(l)}\in\mathbb{R}^{1\times d_{l}}\) are the weights and the biases, respectively; and \(\phi\) is the ReLU activation function. The graph convolutional operator \(A\) is a symmetric normalization of the graph adjacency matrix with self-loops added (Kipf and Welling, 2017).

For ease of exposition, it will be useful to rewrite the matrix notation (1) in element-sums and products. To this end, for a node \(x\), let \(z_{i}^{(l)}(x)\) and \(x_{i}^{(l)}(x)\) denote the pre- and post-activation value at the \(i\)-th coordinate in the \(l\)-th layer, respectively. Particularly, in an \(L\)-layer GCN, \(x^{(0)}(x)\) is the input feature vector and \(z^{(L)}(x)\) is the output vector. The layer architecture of GCN reads

\[y_{i}^{(l)}(x)=\sum_{j=1}^{d_{l-1}}W_{ji}^{(l)}x_{j}^{(l-1)}(x), \quad z_{i}^{(l)}(x)=b_{i}^{(l)}+\sum_{v\in\mathcal{V}}A_{vv}y_{i}^{(l)}(v), \quad x_{i}^{(l)}(x)=\phi(z_{i}^{(l)}(x)). \tag{2}\]

### Limit in the Width

The following theorem states that when the weights and biases in each layer are iid zero-mean Gaussians, in the limit on the layer width, the GCN output \(z^{(L)}(x)\) is a multi-output GP over the index \(x\).

**Theorem 1**.: _Assume \(d_{1},\ldots,d_{L-1}\) to be infinite in succession and let the bias and weight terms be independent with distributions_

\[b_{i}^{(l)}\sim\mathcal{N}(0,\sigma_{b}^{2}),\quad W_{ij}^{(l)} \sim\mathcal{N}(0,\sigma_{w}^{2}/d_{l-1}),\quad l=1,\ldots,L.\]

_Then, for each \(i\), the collection \(\{z_{i}^{(l)}(x)\}\) over all graph nodes \(x\) follows the normal distribution \(\mathcal{N}(0,K^{(l)})\), where the covariance matrix \(K^{(l)}\) can be computed recursively by_

\[C^{(l)}=\mathrm{E}_{z_{i}^{(l)}\sim\mathcal{N}(0,K^{(l)})}[\phi( z_{i}^{(l)})\phi(z_{i}^{(l)})^{T}], l=1,\ldots,L, \tag{3}\] \[K^{(l+1)}=\sigma_{b}^{2}\mathbbm{1}_{N\times N}+\sigma_{w}^{2}AC ^{(l)}A^{T}, l=0,\ldots,L-1. \tag{4}\]

All proofs of this paper are given in the appendix. Note that different from a usual GP, which is a random function defined over a connected region of the Euclidean space, here \(z^{(L)}\) is defined over a discrete set of graph nodes. In the usual use of a graph in machine learning, this set is finite, such that the function distribution degenerates to a multivariate distribution. In semi-supervised learning, the dimension of the distribution, \(N\), is fixed when one conducts transductive learning; but it will vary in the inductive setting because the graph will have new nodes and edges. One special care of a graph-based GP over a usual GP is that the covariance matrix will need to be recomputed from scratch whenever the graph alters.

Theorem 1 leaves out the base definition \(C^{(0)}\), whose entry denotes the covariance between two input nodes. The traditional literature uses the inner product \(C^{(0)}(x,x^{\prime})=\frac{x\cdot x^{\prime}}{d_{0}}\)(Lee et al., 2018),
but nothing prevents us using any positive-definite kernel alternatively.1 For example, we could use the squared exponential kernel \(C^{(0)}(x,x^{\prime})=\exp\left(-\frac{1}{2}\sum_{j=1}^{d_{0}}\left(\frac{x_{j}-x _{j}^{\prime}}{\ell_{j}}\right)^{2}\right)\). Such flexibility in essence performs an implicit feature transformation as preprocessing.

Footnote 1: Here, we abuse the notation and use \(x\) in place of \(x^{(0)}(x)\) in the inner product.

### Universality

A covariance kernel is positive definite; hence, the Moore-Aronszajn theorem (Aronszajn, 1950) suggests that it defines a unique Hilbert space for which it is a reproducing kernel. If this space is dense, then the kernel is called _universal_. One can verify universality by checking if the kernel matrix is positive definite for any set of distinct points.2 For the case of graphs, it suffices to verify if the covariance matrix for all nodes is positive definite.

Footnote 2: Note the conventional confusion in terminology between functions and matrices: a kernel function is positive definite (resp. strictly positive definite) if the corresponding kernel matrix is positive semi-definite (resp. positive definite) for any collection of distinct points.

We do this job for the ReLU activation function. It is known that the kernel \(\mathrm{E}_{w\sim\mathcal{N}(0,I_{d})}[\phi(w\cdot x)\phi(w\cdot x^{\prime})]\) admits a closed-form expression as a function of the angle between \(x\) and \(x^{\prime}\), hence named the _arc-cosine_ kernel (Cho & Saul, 2009). We first establish the following lemma that states that the kernel is universal over a half-space.

**Lemma 2**.: _The arc-cosine kernel is universal on the upper-hemisphere \(S=\left\{x\in\mathbb{R}^{d}:\|x\|_{2}=1,x_{1}>0\right\}\) for all \(d\geq 2\)._

It is also known that the expectation in (3) is proportional to the arc-cosine kernel up to a factor \(\sqrt{K^{(l)}(x,x)K^{(l)}(x^{\prime},x^{\prime})}\)(Lee et al., 2018). Therefore, we iteratively work on the post-activation covariance (3) and the pre-activation covariance (4) and show that the covariance kernel resulting from the limiting GCN is universal, for any GCN with three or more layers.

**Theorem 3**.: _Assume \(A\) is irreducible and non-negative and \(C^{(0)}\) does not contain two linearly dependent rows. Then, \(K^{(0)}\) is positive definite for all \(l\geq 3\)._

### Limit in the Depth

The depth of a neural network exhibits interesting behaviors. Deep learning tends to favor deep networks because of their empirically outstanding performance, exemplified by generations of convolutional networks for the ImageNet classification (Krizhevsky et al., 2012; Wortsman et al., 2022); while graph neural networks are instrumental to be shallow because of the over-smoothing and over-squashing properties (Li et al., 2018; Topping et al., 2022). For multi-layer perceptrons (networks with fully connected layers), several previous works have noted that the recurrence relation of the covariance kernel across layers leads to convergence to a fixed-point kernel, when the depth \(L\rightarrow\infty\) (see, e.g., Lee et al. (2018); in Appendix B.5, we elaborate this limit). In what follows, we offer the parallel analysis for GCN.

**Theorem 4**.: _Assume \(A\) is symmetric, irreducible, aperiodic, and non-negative with Perron-Frobenius eigenvalue \(\lambda>0\). The following results hold as \(l\rightarrow\infty\)._

1. _When_ \(\sigma_{b}^{2}=0\)_,_ \(\rho_{\min}(K^{(l)})\nearrow 1\)_, where_ \(\rho_{\min}\) _denotes the minimum correlation between any two nodes_ \(x\) _and_ \(x^{\prime}\)_._
2. _When_ \(\sigma_{w}^{2}<2/\lambda^{2}\)_, a subsequence of_ \(K^{(l)}\) _converges to some matrix._
3. _When_ \(\sigma_{w}^{2}>2/\lambda^{2}\)_, let_ \(c_{l}=(\sigma_{w}^{2}\lambda^{2}/2)^{l}\)_; then,_ \(K^{(l)}/c_{l}\to vv^{T}\) _where_ \(v\) _is an eigenvector corresponding to_ \(\lambda\)_._

A few remarks follow. The first case implies that the correlation matrix converges monotonously to a matrix of all ones. As a consequence, up to some scaling \(c_{l}^{\prime}\) that may depend on \(l\), the scaled covariance matrix \(K^{(l)}/c_{l}^{\prime}\) converges to a rank-1 matrix. The third case shares a similar result, with the limit explicitly spelled out, but note that the eigenvector \(v\) may not be normalized. The second case is challenging to analyze. According to empirical verification, we speculate a stronger result--convergence of \(K^{(l)}\) to a unique fixed point--may hold.
## 4 Scalable Computation through Low-Rank Approximation

The computation of the covariance matrix \(K^{(L)}\) through recursion (3)-(4) is the main computational bottleneck for GP posterior inference. We start the exposition with the mean prediction. We compute the posterior mean \(\widehat{y}_{*}=K^{(L)}_{sb}(K^{(L)}_{bb}+\epsilon I)^{-1}y_{b}\), where the subscripts \(b\) and \(*\) denote the training set and the prediction set, respectively; and \(\epsilon\), called the _nugget_, is the noise variance of the training data. Let there be \(N_{b}\) training nodes and \(N_{*}\) prediction nodes. It is tempting to compute only the \((N_{b}+N_{*})\times N_{b}\) submatrix of \(K^{(L)}\) for the task, but the recursion (4) requires the full \(C^{(L-1)}\) at the presence of \(A\), and hence all the full \(C^{(l)}\)'s and \(K^{(l)}\)'s.

To reduce the computational costs, we resort to a low-rank approximation of \(C^{(l)}\), from which we easily see that \(K^{(l+1)}\) is also low-rank. Before deriving the approximation recursion, we note (again) that for the ReLU activation \(\phi\), \(C^{(l)}\) in (3) is the arc-cosine kernel with a closed-form expression:

\[C^{(l)}_{xx^{\prime}}=\frac{1}{2\pi}\sqrt{K^{(l)}_{xx}K^{(l)}_{x^{\prime}x^{ \prime}}}\left(\sin\theta^{(l)}_{xx^{\prime}}+(\pi-\theta^{(l)}_{xx^{\prime}} )\cos\theta^{(l)}_{xx^{\prime}}\right)\quad\text{where}\quad\theta^{(l)}_{ xx^{\prime}}=\arccos\left(\frac{K^{(l)}_{xx^{\prime}}}{\sqrt{K^{(l)}_{xx}K^{(l)}_{x^{ \prime}x^{\prime}}}}\right). \tag{5}\]

Hence, the main idea is that starting with a low-rank approximation of \(K^{(l)}\), compute an approximation of \(C^{(l)}\) by using (5), and then obtain an approximation of \(K^{(l+1)}\) based on (4); then, repeat.

To derive the approximation, we use the subscript \(a\) to denote a set of landmark nodes with cardinality \(N_{a}\). The Nystrom approximation (Drineas & Mahoney, 2005) of \(K^{(0)}\) is \(K^{(0)}_{:a}(K^{(0)}_{aa})^{-1}K^{(0)}_{a:}\), where the subscript : denotes retaining all rows/columns. We rewrite this approximation in the Cholesky style as \(Q^{(0)}Q^{(0)}{}^{T}\), where \(Q^{(0)}=K^{(0)}_{:a}(K^{(0)}_{aa})^{-\frac{1}{2}}\) has size \(N\times N_{a}\). Proceed with induction. Let \(K^{(l)}\) be approximated by \(\widehat{K}^{(l)}=Q^{(l)}Q^{(l)}{}^{T}\), where \(Q^{(l)}\) has size \(N\times(N_{a}+1)\). We apply (5) to compute an approximation to \(C^{(l)}_{:a}\), namely \(\widehat{C}^{(l)}_{:a}\), by using \(\widehat{K}^{(l)}_{:a}\). Then, (4) leads to a Cholesky style approximation of \(K^{(l+1)}\):

\[\widehat{K}^{(l+1)}=\sigma_{b}^{2}\mathbf{1}_{N\times N}+\sigma_{w}^{2}A \widehat{C}^{(l)}A^{T}\equiv Q^{(l+1)}Q^{(l+1)}{}^{T},\]

where \(Q^{(l+1)}=\left[\sigma_{w}A\widehat{C}^{(l)}_{:a}(\widehat{C}^{(l)}_{aa})^{- \frac{1}{2}}\quad\sigma_{b}\mathbf{1}_{N\times 1}\right]\). Clearly, \(Q^{(l+1)}\) has size \(N\times(N_{a}+1)\), completing the induction.

In summary, \(K^{(L)}\) is approximated by a rank-\((N_{a}+1)\) matrix \(\widehat{K}^{(L)}=Q^{(L)}{Q^{(L)}}^{T}\). The computation of \(Q^{(L)}\) is summarized in Algorithm 1. Once it is formed, the posterior mean is computed as

\[\widehat{y}_{*}\approx\widehat{K}^{(L)}_{*b}(\widehat{K}^{(L)}_{bb}+\epsilon I )^{-1}y_{b}=Q^{(L)}_{*:}\left({Q^{(L)}_{b:}}^{T}Q^{(L)}_{b:}+\epsilon I\right) ^{-1}{Q^{(L)}_{b:}}^{T}y_{b}, \tag{6}\]

where note that the matrix to be inverted has size \((N_{a}+1)\times(N_{a}+1)\), which is assumed to be significantly smaller than \(N_{b}\times N_{b}\). Similarly, the posterior variance is

\[\widehat{K}^{(L)}_{**}-\widehat{K}^{(L)}_{*b}(\widehat{K}^{(L)}_{bb}+\epsilon I )^{-1}\widehat{K}^{(L)}_{b*}=\epsilon Q^{(L)}_{*:}\left({Q^{(L)}_{b:}}^{T}Q^{ (L)}_{b:}+\epsilon I\right)^{-1}{Q^{(L)}_{*:}}^{T}. \tag{7}\]

The computational costs of \(Q^{(L)}\) and the posterior inference (6)-(7) are summarized in Table 1

\begin{table}
\begin{tabular}{l c c} \hline  & Time \(O(\cdot)\) & Storage \(O(\cdot)\) \\ \hline Computation of \(Q^{(L)}\) & \(LMN_{a}+LNN_{a}^{2}+LN_{a}^{3}\) & \(NN_{a}\) \\ Posterior mean (6) & \(N_{*}N_{a}+N_{b}N_{a}^{2}+N_{a}^{3}\) & \((N_{b}+N_{*})N_{a}\) \\ Posterior variance (7) & \(N_{*}N_{a}+N_{b}N_{a}^{2}+N_{a}^{3}\) & \((N_{b}+N_{*})N_{a}\) \\ \hline \end{tabular}
\end{table}
Table 1: Computational costs. \(M\): number of edges; \(N\): number of nodes; \(N_{b}\): number of training nodes; \(N_{*}\): number of prediction nodes; \(N_{a}\): number of landmark nodes; \(L\): number of layers. Assume \(N_{b}\gg N_{a}\). For posterior variance, assume only the diagonal is needed.
```
0:\(Q^{(0)}\) such that \(K^{(0)}\approx{Q^{(0)}}{Q^{(0)}}^{T}\)
1:for\(l=0,\ldots,L-1\)do
2: Compute \(\widehat{K}^{(l)}_{\alpha}=Q^{(l)}_{:}Q^{(0)}_{a:}\)
3: Compute \(\widehat{C}^{(l)}_{\alpha}\) by (5), where \(C^{(l)}\) (resp. \(K^{(l)}\)) entries are replaced by \(\widehat{C}^{(l)}\) (resp. \(\widehat{K}^{(l)}\)) entries
4: Compute \(Q^{(l+1)}=\left[\sigma_{w}\widehat{A}\widehat{C}^{(l)}_{\alpha a}(\widehat{C} ^{(l)}_{aa})^{-\frac{1}{2}}\quad\sigma_{b}\mathbf{1}_{N\times 1}\right]\)
5:endfor
```

**Algorithm 1** Computing \(K^{(L)}\approx\widehat{K}^{(L)}={Q^{(L)}}{Q^{(L)}}^{T}\)

## 5 Composing Graph Neural Network-Inspired Kernels

Theorem 1, together with its proof, suggests that the covariance matrix of the limiting GP can be computed in a composable manner. Moreover, the derivation of Algorithm 1 indicates that the low-rank approximation of the covariance matrix can be similarly composed. Altogether, such a nice property allows one to easily derive the corresponding covariance matrix and its approximation for a new GNN architecture, like writing a program and obtaining a transformation of it automatically through operator overloading (Novak et al., 2020): the covariance matrix is a transformation of the GNN and the composition of the former is in exactly the same manner and order as that of the latter. We call the covariance matrices _programmable_.

For example, we write a GCN layer as \(X\gets A\phi(X)W+b\), where for notational simplicity, \(X\) denotes pre-activation rather than post-activation as in the earlier sections. The activation \(\phi\) on \(X\) results in a transformation of the kernel matrix \(K\) into \(g(K)\), defined as:

\[g(K):=C=\mathrm{E}_{z\sim\mathcal{N}(0,K)}[\phi(z)\phi(z)^{T}], \tag{8}\]

due to (3). Moreover, if \(K\) admits a low-rank approximation \(QQ^{T}\), then \(g(K)\) admits a low-rank approximation \(PP^{T}\) where \(P=\mathrm{Chol}(g(K))\) with

\[\mathrm{Chol}(C):=C_{:a}C_{aa}^{-\frac{1}{2}}.\]

The next operation--graph convolution--multiplies \(A\) to the left of the post-activation. Correspondingly, the covariance matrix \(K\) is transformed to \(AKA^{T}\) and the low-rank approximation factor \(Q\) is transformed to \(AQ\). Then, the operation--multiplying the weight matrix \(W\) to the right--will transform \(K\) to \(\sigma_{w}^{2}K\) and \(Q\) to \(\sigma_{w}Q\). Finally, adding the bias \(b\) will transform \(K\) to \(K+\sigma_{b}^{2}\mathbf{1}_{N\times N}\) and \(Q\) to \([Q\quad\sigma_{b}\mathbf{1}_{N\times 1}]\). Altogether, we have obtained the following updates per layer:

\[\text{GCN}:\quad X \gets A\phi(X)W+b\] \[K \leftarrow\sigma_{w}^{2}Ag(K)A^{T}+\sigma_{b}^{2}\mathbf{1}_{N \times N}\] \[Q \leftarrow\left[\sigma_{w}A\,\mathrm{Chol}(g(QQ^{T}))\quad\sigma_{ b}\mathbf{1}_{N\times 1}\right].\]

One may verify the \(K\) update against (3)-(4) and the \(Q\) update against Algorithm 1. Both updates can be automatically derived based on the update of \(X\).

We summarize the building blocks of a GNN and the corresponding kernel/low-rank operations in Table 2. The independent-addition building block is applicable to skip/residual connections. For

\begin{table}
\begin{tabular}{l l l l} \hline Building block & Neural network & Kernel operation & Low-rank operation \\ \hline Input & \(X\gets X^{(0)}\) & \(K\gets C^{(0)}\) & \(Q\leftarrow\mathrm{Chol}(C^{(0)})\) \\ Bias term & \(X\gets X+b\) & \(K\gets K+\sigma_{b}^{2}\mathbf{1}_{N\times N}\) & \(Q\leftarrow[Q\quad\sigma_{b}\mathbf{1}_{N\times 1}]\) \\ Weight term & \(X\gets XW\) & \(K\gets\sigma_{w}^{2}K\) & \(Q\leftarrow\sigma_{w}Q\) \\ Mixed weight term & \(X\gets X(\alpha I+\beta W)\) & \(K\leftarrow(\alpha^{2}+\beta^{2}\sigma_{w}^{2})K\) & \(Q\leftarrow\sqrt{\alpha^{2}+\beta^{2}\sigma_{w}^{2}}Q\) \\ Graph convolution & \(X\gets AX\) & \(K\gets AKA^{T}\) & \(Q\gets AQ\) \\ Activation & \(X\leftarrow\phi(X)\) & \(K\gets g(K)\) & \(Q\leftarrow\mathrm{Chol}(g(QQ^{T}))\) \\ Independent addition & \(X\gets X_{1}+X_{2}\) & \(K\gets K_{1}+K_{2}\) & \(Q\leftarrow[Q_{1}\quad Q_{2}]\) \\ \hline \end{tabular}
\end{table}
Table 2: Neural network building blocks, kernel operations, and the low-rank counterpart.
example, here is the composition for the GCNII layer (Chen et al., 2020) without a bias term, where a skip connection with \(X^{(0)}\) occurs:

\[\text{GCNII}: X\leftarrow\left((1-\alpha)A\phi(X)+\alpha X^{(0)}\right)((1- \beta)I+\beta W)\] \[K\leftarrow\left((1-\alpha)^{2}Ag(K)A^{T}+\alpha^{2}K^{(0)} \right)((1-\beta)^{2}+\beta^{2}\sigma_{w}^{2})\] \[Q\leftarrow\left[(1-\alpha)A\operatorname{Chol}(g(QQ^{T})) \quad\alpha Q^{(0)}\right]\sqrt{(1-\beta)^{2}+\beta^{2}\sigma_{w}^{2}}.\]

For another example of the composability, we consider the popular GIN layer (Xu et al., 2019), which we assume uses a 2-layer MLP after the neighborhood aggregation:

\[\text{GIN}: X\leftarrow\phi(A\phi(X)W+b)W^{\prime}+b^{\prime}\] \[K\leftarrow\sigma_{w}^{2}g(B)+\sigma_{b}^{2}\mathbf{1}_{N\times N }\quad\text{where}\quad B=\sigma_{w}^{2}Ag(K)A^{T}+\sigma_{b}^{2}\mathbf{1}_{ N\times N}\] \[Q\leftarrow\left[\sigma_{w}\operatorname{Chol}(g(PP^{T}))\quad \sigma_{b}\mathbf{1}_{N\times 1}\right]\quad\text{where}\quad P=\left[\sigma_{w}A \operatorname{Chol}(g(QQ^{T}))\quad\sigma_{b}\mathbf{1}_{N\times 1}\right].\]

Additionally, the updates for a GraphSAGE layer (Hamilton et al., 2017) are given in Appendix C.

## 6 Experiments

In this section, we conduct a comprehensive set of experiments to evaluate the performance of the GP kernels derived by taking limits on the layer width of GCN and other GNNs. We demonstrate that these GPs are comparable with GNNs in prediction performance, while being significantly faster to compute. We also show that the low-rank version scales favorably, suitable for practical use.

**Datasets.** The experiments are conducted on several benchmark datasets of varying sizes, covering both classification and regression. They include predicting the topic of scientific papers organized in a citation network (Cora, Citeseer, PubMed, and ArXiv); predicting the community of online posts based on user comments (Reddit), and predicting the average daily traffic of Wikipedia pages using hyperlinks among them (Chameleon, Squirrel, and Crocodile). Details of the datasets (including sources and preprocessing) are given in Appendix D.

**Experiment environment, training details, and hyperparameters** are given in Appendix D.

**Prediction Performance: GCN-based comparison.** We first conduct the semi-supervised learning tasks on all datasets by using GCN and GPs with different kernels. These kernels include the one equivalent to the limiting GCN (GCNGP), a usual squared-exponential kernel (RBF), and the GGP kernel proposed by Ng et al. (2018).3 Each of these kernels has a low-rank version (suffixed with -X). RBF-X and GGP-X4 use the Nystrom approximation, consistent with GCNGP-X.

Footnote 3: We apply only the kernel but not the likelihood nor the variational inference used in Ng et al. (2018), for reasons given in Appendix D.

Footnote 4: GGP-X in our notation is the Nyström approximation of the GGP kernel, different from a technique under the same name in Ng et al. (2018), which uses additionally the validation set to compute the prediction loss.

GPs are by nature suitable for regression. For classification tasks, we use the one-hot representation of labels to set up a multi-output regression. Then, we take the coordinate with the largest output as the class prediction. Such an ad hoc treatment is widely used in the literature, as other more principled approaches (such as using the Laplace approximation on the non-Gaussian posterior) are too time-consuming for large datasets, meanwhile producing no noticeable gain in accuracy.

Table 3 summarizes the accuracy for classification and the coefficient of determination, \(R^{2}\), for regression. Whenever randomness is involved, the performance is reported as an average over five runs. The results of the two tasks show different patterns. For classification, GCNGP-(X) is sligtly better than GCN and GGP-(X), while RBF-(X) is significantly worse than all others; moreover, the low-rank version is outperformed by using the full kernel matrix. On the other hand, for regression, GCNGP-(X) significantly outperforms GCN, RBF-(X), and GGP-(X); and the low-rank version becomes better. The less competitive performance of RBF-(X) is expected, as it does not leverage the graph inductive bias. It is attractive that GCNGP-(X) is competitive with GCN.

**Prediction Performance: Comparison with other GNNs.** In addition to GCN, we conduct experiments with several popularly used GNN architectures (GCNII, GIN, and GraphSAGE) and GPs
with the corresponding kernels. We test with the three largest datasets: PubMed, ArXiv, and Reddit, for the latter two of which a low-rank version of the GPs is used for computational feasibility.

Table 4 summarizes the results. The observations on other GNNs extend similarly those on the GCN. In particular, on PubMed the GPs noticeably improve over the corresponding GNNs, while on ArXiv and Reddit the two families perform rather similarly. An exception is GIN for ArXiv, which significantly underperforms the GP counterpart, as well as other GNNs. It may improve with an extensive hyperparameter tuning.

**Running time.** We compare the running time of the methods covered by Table 3. Different from usual neural networks, the training and inference of GNNs do not decouple in full-batch training. Moreover, there is not a universally agreed split between the training and the inference steps in GPs. Hence, we compare the total time for each method.

Figure 1 plots the timing results, normalized against the GCN time for ease of comparison. It suggests that GCNGP(-X) is generally faster than GCN. Note that the vertical axis is in the logarithmic scale. Hence, for some of the datasets, the speedup is even one to two orders of magnitude.

**Scalability.** For graphs, especially under the semi-supervised learning setting, the computational cost of a GP is much more complex than that of a usual one (which can be simply described as "cubic in the training set size"). One sees in Table 1 the many factors that determine the cost of our graph-based low-rank kernels. To explore the practicality of the proposed method, we use the timings gathered for Figure 1 to obtain an empirical scaling with respect to the graph size, \(M+N\).

Figure 2 fits the running times, plotted in the log-log scale, by using a straight line. We see that for neither GCN nor GCNGP(-X), the actual running time closely follows a polynomial complexity. However, interestingly, the least-squares fittings all lead to a slope approximately 1, which agrees with a linear cost. Theoretically, only GCNGP-X and GCN are approximately linear with respect to \(M+N\), while GCNGP is cubic.

**Analysis on the depth.** The performance of GCN deteriorates with more layers, known as the oversmoothing phenomenon. Adding residual/skip connections mitigates the problem, such as in GCNII. A natural question asks if the corresponding GP kernels behave similarly.

Figure 3 shows that the trends of GCN and GCNII are indeed as expected. Interestingly, their GP counterparts both remain stable for depth \(L\) as large as 12. Our depth analysis (Theorem 4) suggests that in the limit, the GPs may perform less well because the kernel matrix may degenerate to rank 1. This empirical result indicates that the drop in performance may have not started yet.

\begin{table}
\begin{tabular}{c|c c|c c|c c} \hline \hline  & \multicolumn{2}{c|}{PubMed} & \multicolumn{2}{c|}{ArXiv} & \multicolumn{2}{c}{Reddit} \\ architecture & GNN & GNNGP & GNN & GNNGP-X & GNN & GNNGP-X \\ \hline GCN & 0.7649\(\pm\)0.0058 & 0.7960 & 0.6989\(\pm\)0.0016 & 0.7011\(\pm\)0.0011 & 0.9330\(\pm\)0.0006 & 0.9465\(\pm\)0.0003 \\ GCNII & 0.7558\(\pm\)0.0096 & 0.7840 & 0.7008\(\pm\)0.0021 & 0.6955\(\pm\)0.0011 & 0.9482\(\pm\)0.0007 & 0.9500\(\pm\)0.0003 \\ GIN & 0.7406\(\pm\)0.0112 & 0.7690 & 0.6340\(\pm\)0.0066 & 0.6652\(\pm\)0.0012 & 0.9398\(\pm\)0.0016 & 0.9428\(\pm\)0.0005 \\ GraphSAGE & 0.7535\(\pm\)0.0047 & 0.7900 & 0.6984\(\pm\)0.0021 & 0.6962\(\pm\)0.0007 & 0.9628\(\pm\)0.0007 & 0.9539\(\pm\)0.0003 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Performance comparison (Micro-F1) between GNNs and the corresponding GP kernels.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline  & GCN & GCNGP & GCNGP-X & RBF & RBF-X & GGP & GGP-X \\ \hline Cora & 0.8183\(\pm\)0.0055 & **0.8280** & 0.7980 & 0.5860 & 0.5850 & 0.7850 & 0.7410 \\ Citeseer & 0.6941\(\pm\)0.0079 & **0.7090** & 0.7080 & 0.6120 & 0.6090 & 0.7060 & 0.6470 \\ PubMed & 0.7649\(\pm\)0.0058 & **0.7960** & 0.7810 & 0.7360 & 0.7340 & 0.7820 & 0.7380 \\ ArXiv & 0.6990\(\pm\)0.0014 & OOM & **0.7011\(\pm\)0.0011** & OOM & 0.5382 & OOM & 0.6527 \\ Reddit & 0.9330\(\pm\)0.0006 & OOM & **0.9465\(\pm\)0.0003** & OOM & 0.5920 & OOM & 0.9058 \\ \hline Chameleon & 0.5690\(\pm\)0.0376 & 0.6720 & **0.6852** & 0.5554 & 0.5613 & 0.5280 & 0.5311 \\ Squirrel & 0.4243\(\pm\)0.0393 & 0.4926 & **0.4998** & 0.3187 & 0.3185 & 0.2440 & 0.2251 \\ Crocodile & 0.6976\(\pm\)0.0323 & 0.8002 & **0.8013** & 0.6643 & 0.6710 & 0.6952 & 0.6810 \\ \hline \hline \end{tabular}
\end{table}
Table 3: Performance of GCNGP, in comparison with GCN and typical GP kernels. The Micro-F1 score is reported for classification tasks and \(R^{2}\) is reported for regression tasks.
**Analysis on the landmark set.** The number of landmark nodes, \(N_{a}\), controls the approximation quality of the low-rank kernels and hence the prediciton accuracy. On the other hand, the computational costs summarized in Table 1 indicate a dependency on \(N_{a}\) as high as the third power. It is crucial to develop an empirical understanding of the accuracy-time trade-off it incurs.

Figure 4 clearly shows that as \(N_{a}\) becomes larger, the running time increase is not linear, while the increase of accuracy diminishes as the landmark set approaches the training set. It is remarkable that using only 1/800 of the training set as landmark nodes already achieves an accuracy surpassing that of GCN, by using time that is only a tiny fraction of the time otherwise needed to gain an additional 1% increase in the accuracy.

## 7 Conclusions

We have presented a GP approach for semi-supervised learning on graph-structured data, where the covariance kernel incorporates a graph inductive bias by exploiting the relationship between a GP and a GNN with infinitely wide layers. Similar to other neural networks priorly investigated, one can work out the equivalent GP (in particular, the covariance kernel) for GCN; and inspired by this equivalence, we formulate a procedure to compose covariance kernels corresponding to many other members of the GNN family. Moreover, every building block in the procedure has a low-rank counterpart, which allows for building a low-rank approximation of the covariance matrix that facilitates scalable posterior inference. We demonstrate the effectiveness of the derived kernels used for semi-supervised learning and show their advantages in computation time over GNNs.
## Acknowledgments

Mihai Anitescu was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research (ASCR) under Contract DE-AC02-06CH11347. Jie Chen acknowledges supports from the MIT-IBM Watson AI Lab.
