LazyGNN: Large-Scale Graph Neural Networks via Lazy Propagation
Recent works have demonstrated the benefits of capturing long-distance
dependency in graphs by deeper graph neural networks (GNNs). But deeper GNNs
suffer from the long-lasting scalability challenge due to the neighborhood
explosion problem in large-scale graphs. In this work, we propose to capture
long-distance dependency in graphs by shallower models instead of deeper
models, which leads to a much more efficient model, LazyGNN, for graph
representation learning. Moreover, we demonstrate that LazyGNN is compatible
with existing scalable approaches (such as sampling methods) for further
accelerations through the development of mini-batch LazyGNN. Comprehensive
experiments demonstrate its superior prediction performance and scalability on
large-scale benchmarks. The implementation of LazyGNN is available at
https://github.com/RXPHD/Lazy_GNN.
Rui Xue, Haoyu Han, MohamadAli Torkamani, Jian Pei, Xiaorui Liu
# LazyGNN: Large-Scale Graph Neural Networks via Lazy Propagation

###### Abstract

Recent works have demonstrated the benefits of capturing long-distance dependency in graphs by deeper graph neural networks (GNNs). But deeper GNNs suffer from the long-lasting scalability challenge due to the neighborhood explosion problem in large-scale graphs. In this work, we propose to capture long-distance dependency in graphs by shallower models instead of deeper models, which leads to a much more efficient model, LazyGNN, for graph representation learning. Moreover, we demonstrate that LazyGNN is compatible with existing scalable approaches (such as sampling methods) for further accelerations through the development of mini-batch LazyGNN. Comprehensive experiments demonstrate its superior prediction performance and scalability on large-scale benchmarks. LazyGNN also achieves state-of-art performance on the OGB leaderboard.

Machine Learning, ICML

## 1 Introduction

Graph neural networks (GNNs) have been widely used for representation learning on graph-structured data (Hamilton, 2020; Ma and Tang, 2021), and they achieve promising state-of-the-art performance on various general graph learning tasks, such as node classification, link prediction, and graph classification (Kipf and Welling, 2016; Gasteiger et al., 2019; Velickovic et al., 2017; Wu et al., 2019) as well as a variety of important applications, such as recommendation systems, social network analysis, and transportation prediction. In particular, recent research in deeper GNNs has generally revealed the performance gains from capturing long-distance relations in graphs by stacking more graph convolution layers or unrolling various fixed point iterations (Gasteiger et al., 2018; Gu et al., 2020; Liu et al., 2020; Chen et al., 2020; Li et al., 2021; Ma et al., 2020; Pan et al., 2020; Zhu et al., 2021; Chen et al., 2020). However, the recursive feature propagations in deeper GNNs lead to the well-known neighborhood explosion problem since the number of neighbors grows exponentially with the number of feature propagation layers (Hamilton et al., 2017; Chen et al., 2018). This causes tremendous scalability challenges for data sampling, computation, memory, parallelism, and end-to-end training when employing GNNs on large-scale graphs. It greatly limits GNNs' broad applications in large-scale industry-level applications due to limited computation and memory resources (Ying et al., 2018; Shao et al., 2022).

A large body of existing research focuses on improving the scalability and efficiency of large-scale GNNs using various innovative designs, such as sampling methods, pre-computing or post-computing methods, and distributed methods. Although these approaches mitigate the neighborhood explosion problem, they still face various limitations when they are applied to deeper GNNs. For instance, sampling approaches (Hamilton et al., 2017; Chen et al., 2018; Zeng et al., 2020; Zou et al., 2019; Fey et al., 2021; Yu et al., 2022) usually incur large approximation errors and suffer from performance degradation as demonstrated in large-scale OGB benchmarks or require large additional memory or storage space to store activation values in hidden layers. Pre-computing or post-computing methods (Wu et al., 2019; Rossi et al., 2020; Sun et al., 2021; Zhang et al., 2022; Bojchevski et al., 2020; Zhu, 2005; Huang et al., 2020) lose the benefits of end-to-end training and usually suffer from performance loss. Distributed methods (Chiang et al., 2019; Chai et al., 2022; Shao et al., 2022) distribute large graphs to multiple servers for parallel training, but they either lose the inter-server edges or suffer from expensive feature communication between servers.

In this work, we take a substantially different and novel perspective and propose to capture long-distance dependency in graphs by shallower GNNs instead of deeper GNNs. The key intuition comes from the fact that node information has been propagated over the graph many times during the training process so we may only need to propagate information lazily by reusing the propagated information over the training iterations. This intuition leads to the proposed LazyGNN that solves the inherent neighborhood explosion
problem by significantly reducing the number of aggregation layers while still capturing long-distance dependency in graphs through lazy propagation. Multiple technical challenges such as the risk of over-smoothing, additional variation due to feature dropout, and back-propagation through historical computation graphs are carefully dealt with by innovative designs in LazyGNN. Moreover, since LazyGNN is a shallow model, it naturally tackles the limitations that existing scalable approaches suffer from when handling large-scale and deeper GNNs. Therefore, various scalable approaches can be used in LazyGNN for further acceleration. We demonstrate this by developing a highly scalable and efficient mini-batch LazyGNN based on sampling methods. To the best of our knowledge, LazyGNN is the first GNN model being versatilely friendly to data sampling, computation, memory space, parallelism, and end-to-end training but still captures long-distance dependency in graphs. Our contributions can be summarized as follows:

* We reveal a novel perspective to solve the neighborhood explosion problem by exploiting the computation redundancy in training GNNs;
* A novel shallow model, LazyGNN, is proposed to capture long-distance dependency in graphs for end-to-end graph representation learning through lazy forward and backward propagations;
* We demonstrate that existing scalable approaches can be compatible with LazyGNN by introducing a highly efficient and scalable mini-batch LazyGNN to handle large-scale graphs based on sampling methods.
* Comprehensive experiments demonstrate superior prediction performance and efficiency on large-scale graph datasets. Notably, the proposed LazyGNN also achieves the best performance on OGB leaderboard.

## 2 Preliminaries

**Notations.** A graph is represented by \(\mathcal{G}=(\mathcal{V},\mathcal{E})\) where \(\mathcal{V}=\{v_{1},\ldots,v_{N}\}\) is the set of nodes and \(\mathcal{E}=\{e_{1},\ldots,e_{M}\}\) is the set of edges. Suppose that each node is associated with a \(d\)-dimensional feature vector, and the original features for all nodes are denoted as \(\mathbf{X}_{\text{fea}}\in\mathbb{R}^{N\times d}\). The graph structure of \(\mathcal{G}\) can be represented by an adjacency matrix \(\mathbf{A}\in\mathbb{R}^{N\times N}\), where \(\mathbf{A}_{ij}>0\) when there exists an edge between node \(v_{i}\) and \(v_{j}\), otherwise \(\mathbf{A}_{i,j}=0\). The symmetrically normalized graph Laplacian matrix is defined as \(\tilde{\mathbf{L}}=\mathbf{I}-\tilde{\mathbf{A}}\) with \(\tilde{\mathbf{A}}=\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}\) where \(\mathbf{D}\) is the degree matrix. Next, we briefly introduce the graph signal denoising and fixed point iteration perspective of GNNs that provide a better understanding of the computation trajectory and long-distance dependency in graphs. Finally, we provide a preliminary study to reveal the computation redundancy in GNNs.

### GNNs as Graph Signal Denoising

It is recently established that many popular GNNs can be uniformly understood as solving graph signal denoising problems with various diffusion properties and that the long-distance dependency can be well captured by unrolling various fixed point iterations (Ma et al., 2020; Pan et al., 2020; Zhu et al., 2021; Chen et al., 2020; Gu et al., 2020). For instance, the message passing in GCN (Kipf and Welling, 2016), \(\mathbf{X}_{\text{out}}=\tilde{\mathbf{A}}\mathbf{X}_{\text{in}}\), can be considered as one gradient descent iteration for minimizing \(\operatorname{tr}(\mathbf{X}^{\top}(\mathbf{I}-\tilde{\mathbf{A}})\mathbf{X})\) with the initialization \(\mathbf{X}_{0}=\mathbf{X}_{\text{in}}\). The message passing scheme in APPNP (Klicpera et al., 2018) follows the aggregation \(\mathbf{X}_{l+1}=(1-\alpha)\tilde{\mathbf{A}}\mathbf{X}_{l}+\alpha\mathbf{X}_ {\text{in}}\) that iteratively minimizes \(\|\mathbf{X}-\mathbf{X}_{\text{in}}\|_{F}^{2}+(1/\alpha-1)\,\operatorname{tr }(\mathbf{X}^{\top}(\mathbf{I}-\tilde{\mathbf{A}})\mathbf{X})\) with the initialization \(\mathbf{X}_{0}=\mathbf{X}_{\text{in}}\) where \(l\) is the index of layers. Implicit GNN (Gu et al., 2020) adopts projected gradient descent to solve the fixed point problem \(\mathbf{X}=\phi(\mathbf{W}\mathbf{X}\tilde{\mathbf{A}}+\mathbf{B})\). Please refer to the reference (Ma et al., 2020; Pan et al., 2020; Zhu et al., 2021; Chen et al., 2020) for such understanding of many other popular GNN models. Moreover, a large number of advanced GNN models have been inspired from this perspective (Chen et al., 2022; Liu et al., 2021; Yang et al., 2021; Jia and Benson, 2022; Jiang et al., 2022).

### Computation Redundancy in Training GNNs

In the training process of GNNs, the model parameters are updated following gradient descent style algorithms such as Adam (Kingma and Ba, 2014). Therefore, the model as well as hidden features in GNNs changes smoothly, especially in the late stages when both the learning rate and gradient norm diminish. This intuition motivates us to investigate the computation redundancy in GNNs. Specifically, we measure the relative hidden features changes, i.e., \(\|\mathbf{X}^{k+1}-\mathbf{X}^{k}\|_{F}/\|\mathbf{X}^{k}\|_{F}\), over the training iterations on Cora dataset using APPNP. We show two cases in Figure 1: (1) when there is no dropout, the hidden features barely change as the training goes; (2) if we use dropout, it will incur additional variation in the hidden features due to the randomness in dropout layers. Both cases demonstrate that the activation values in hidden layers of GNNs indeed change slowly, indicating the existence of _computation redundancy_: the computation in successive training iterations is highly similar. This observation not only validates the rationality of historical embedding used in existing works such as VRGNN (Chen et al., 2017) and GNNAutoScale (Fey et al., 2021) but also motivates the lazy propagation in this work.

## 3 Lazy Graph Neural Networks

In this section, we propose a novel shallow LazyGNN that uses a few aggregation layers to capture long-distance dependency in graphs through lazy forward and backward propagations. Then a mini-batch LazyGNN is introduced to
handle large-scale graphs with efficient data sampling, computation, and memory usage. Complexity analyses are also provided to illustrate the superior scalability of LazyGNN.

### GNNs with Lazy Propagation

Existing research has demonstrated the benefits of capturing long-distance relations in graphs by stacking more feature aggregation layers or unrolling various fixed point iterations as introduced in Section 1. However, these deeper GNNs are not efficient due to the neighborhood explosion problem. Our preliminary study in Section 2 reveals the computation redundancy in GNNs over the training iterations: the hidden features in GNNs evolve slowly such that the computation of feature aggregations is highly correlated and redundant over training iterations. This observation motivates us to develop a novel shallow LazyGNN that captures long-distance relations in graphs by propagating information lazily with very few message-passing layers.

Without loss of generality, we build LazyGNN on the most common graph signal denoising problem:

\[\min_{\mathbf{X}}\ \|\mathbf{X}-\mathbf{X}_{\text{in}}\|_{F}^{2}+(1/\alpha-1)\ \mathrm{tr}(\mathbf{X}^{\top}(\mathbf{I}-\tilde{\mathbf{A}})\mathbf{X}), \tag{1}\]

where the first term maintains the proximity with node hidden features \(\mathbf{X}_{\text{in}}\) after feature transformation, and the second Laplacian smoothing regularization encodes smoothness assumption on graph representations. We take this as an example to illustrate the main idea of lazy propagation but the idea can be applied to other GNNs inspired by general denoising problems or fixed point iterations with different properties. Next, we illustrate the lazy forward and backward propagations in LazyGNN as well as the technical challenges being solved by innovative designs.

**Forward Computation.** From Eq. (1), we derive the high-order graph diffusion as in APPNP (Klicpera et al., 2018) with \(f\) being the feature transformation:

\[\mathbf{X}_{\text{in}}^{k} =f(\mathbf{X}_{\text{fea}},\Theta^{k}), \tag{2}\] \[\mathbf{X}_{0}^{k} =\mathbf{X}_{\text{in}}^{k},\] (3) \[\mathbf{X}_{l+1}^{k} =(1-\alpha)\tilde{\mathbf{A}}\mathbf{X}_{l}^{k}+\alpha\mathbf{X}_ {\text{in}}^{k},\forall l=0,\dots,L-1 \tag{4}\]

where \(l\) and \(k\) denote the index of layers and training iterations, respectively. The key insight of LazyGNN is that the approximate solution of Eq. (1) (i.e., \(\mathbf{X}_{L}^{k}\)) evolves smoothly since \(\mathbf{X}_{\text{in}}^{k}=f(\mathbf{X}_{\text{fea}},\Theta^{k})\) changes smoothly with model parameters \(\Theta^{k}\). Intuitively, the features have been propagated over the graph many times in previous training iterations so that it suffices to propagate features lazily by reusing existing computation. Formally, we propose LazyGNN to leverage the computation redundancy between training iterations by mixing the diffusion output in iteration \(k-1\) (i.e., \(\mathbf{X}_{L}^{k-1}\)) into the initial embedding of the diffusion process in training iteration \(k\), namely \(\mathbf{X}_{0}^{k}=(1-\beta)\mathbf{X}_{L}^{k-1}+\beta\mathbf{X}_{\text{in}}^ {k}\). This is because, in practice, dropout is commonly used in deep learning to prevent overfitting, which introduces additional variations in the feature embedding as demonstrated in Figure 1. Therefore, we introduce this momentum correction to compensate for such disturbance by mixing current and history features with hyperparameter \(\beta\). To summarize, the forward computation of LazyGNN works as follows:

\[\mathbf{X}_{\text{in}}^{k} =f(\mathbf{X}_{\text{fea}},\Theta^{k}), \tag{5}\] \[\mathbf{X}_{0}^{k} =(1-\beta)\mathbf{X}_{L}^{k-1}+\beta\mathbf{X}_{\text{in}}^{k},\] (6) \[\mathbf{X}_{l+1}^{k} =(1-\alpha)\tilde{\mathbf{A}}\mathbf{X}_{l}^{k}+\alpha\mathbf{X} _{\text{in}}^{k},\forall l=0,\dots,L-1 \tag{7}\]

By lingering the computation over training iterations as shown in Figure 2, the feature aggregation layers in Eq. (7) solve the denoising problem in Eq. (1) with an implicitly large number of steps although there are only \(L\) layers in each training iteration. Therefore, it only requires a few feature aggregation layers (a small \(L\)) to approximate the fixed point solution \(\mathbf{X}_{*}^{k}\). Note that we find \(L=1\) and \(L=2\) work very well in our experiments in Section 4. In the extreme case when the learning rate and dropout rate are \(0\), the forward computation of LazyGNN is equivalent to running forward propagations \(L\times K\) times continuously with \(K\) being the total number of training iterations.

**Remark 1** (Long-distance dependency).: _Through lazy propagation, LazyGNN is able to capture long-distance dependency in graphs with a small number of feature aggregation layers because LazyGNN accumulatively propagates information to distant nodes in the graphs over many training iterations. In contrast, existing works try to capture long-distance dependency in graphs by increasing the number of feature propagation layers, which is less efficient._

Figure 2: LazyGNN with Lazy forward & backward propagation.
**Remark 2** (Over-smoothing).: _In contrast to many GNN models such as GCN or GAT that suffer from the over-moothing issue when more layers are stacked, the accumulation of feature aggregations over training iterations in LazyGNN will not cause the over-smoothing issue because the residual connection \(\mathbf{X}_{\text{in}}^{k}\) in Eq. (7) determines the fixed point and prevents the over-smoothed solution._

**Backward Computation.** While it is feasible to leverage the computation from previous training iterations in the forward computation, it is highly non-trivial to compute the gradient for model update in the backward computation since the computation graphs from previous training iterations have been destroyed and released in the memory. In other words, there is no way to directly compute the backpropagation through the history variables \(\{\mathbf{X}_{L}^{k-1},\mathbf{X}_{L}^{k-2},\dots\}\) in current iterations \(k\). Therefore, we choose to compute the gradient indirectly via the implicit function theorem.

**Theorem 1** (Implicit Gradient).: _Let \(\mathbf{X}_{*}\) be the fixed point solution of function \(g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})\), i.e., \(g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})=\mathbf{0}\). Given the gradient of loss function \(\mathcal{L}(\mathbf{X}_{*},\mathbf{Y})\) with respect to the fixed point \(\mathbf{X}_{*}\), i.e., \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{*}}\), the gradient of loss \(\mathcal{L}\) with respect to feature \(\mathbf{X}_{\text{in}}\) can be computed as:_

\[\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}}=\frac{\partial \mathcal{L}}{\partial\mathbf{X}_{*}}(\mathbf{J}|_{\mathbf{X}_{*}})^{-1}\frac{ \partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial\mathbf{X}_{\text{ in}}} \tag{8}\]

_where \(\mathbf{J}|_{\mathbf{X}_{*}}=\frac{\partial g(\mathbf{X}_{*},\mathbf{X}_{ \text{in}})}{\partial\mathbf{X}_{*}}\) is the Jacobian matrix of \(g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})\) evaluated at \(\mathbf{X}_{*}\)._

Proof.: Take the first-order derivative on both sides of the fixed point equation \(g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})=0\):

\[\frac{\partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial \mathbf{X}_{\text{in}}}+\frac{\partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}} )}{\partial\mathbf{X}_{*}}\frac{d\mathbf{X}_{*}}{d\mathbf{X}_{\text{in}}}= \mathbf{0}\] \[\frac{d\mathbf{X}_{*}}{d\mathbf{X}_{\text{in}}}=-\Big{(}\frac{ \partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial\mathbf{X}_{*}} \Big{)}^{-1}\frac{\partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial \mathbf{X}_{\text{in}}}\]

Using \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}}=\frac{\partial \mathcal{L}}{\partial\mathbf{X}_{*}}\frac{d\mathbf{X}_{*}}{d\mathbf{X}_{\text{ in}}}\), we obtain:

\[\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}}=-\frac{\partial \mathcal{L}}{\partial\mathbf{X}_{*}}\Big{(}\frac{\partial g(\mathbf{X}_{*}, \mathbf{X}_{\text{in}})}{\partial\mathbf{X}_{*}}\Big{)}^{-1}\frac{\partial g (\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial\mathbf{X}_{\text{in}}}.\]

Specifically, for the problem in Eq. (1), we have the fixed point equation:

\[g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})=\mathbf{X}_{*}-\mathbf{X}_{\text{in }}+(\frac{1}{\alpha}-1)(\mathbf{I}-\tilde{\mathbf{A}})\mathbf{X}_{*}=\mathbf{0} \tag{9}\]

which gives \(\frac{\partial g(\mathbf{X}_{*},\mathbf{X}_{\text{in}})}{\partial\mathbf{X}_{ \text{in}}}=-\mathbf{I}\) and \(\mathbf{J}|_{\mathbf{X}_{*}}=\frac{1}{\alpha}(\mathbf{I}-(1-\alpha)\tilde{ \mathbf{A}})\). Therefore, according to Theorem 1, the gradient of loss \(\mathcal{L}\) with respect to \(\mathbf{X}_{\text{in}}\) can be computed as follows:

\[\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}}=\alpha\frac{ \partial\mathcal{L}}{\partial\mathbf{X}_{*}}\Big{(}\mathbf{I}-(1-\alpha)\tilde{ \mathbf{A}}\Big{)}^{-1}. \tag{10}\]

However, it is still infeasible to compute the expensive matrix inversion so we propose to approximate it by the iterative backward gradient propagation:

\[\mathbf{G}_{L}=\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}} \leavevmode\nobreak\ \Big{(}\approx\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{*}}\Big{)} \tag{11}\] \[\mathbf{G}_{l}=(1-\alpha)\tilde{\mathbf{A}}\mathbf{G}_{l+1}+\alpha \frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}}\leavevmode\nobreak\ \forall l=L-1,\dots,0 \tag{12}\]

where \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{*}}\) is approximated by \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}}\), and \(\mathbf{G}_{0}\) provides an approximation for gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}}\). It is clear that the backward computation requires gradient propagation over the graph, which is symmetric to and as expensive as the vanilla forward computation. Similarly, to reduce the number of gradient propagation layers, we propose to propagate the gradient lazily by leveraging the gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}^{k-1}}\) computed in the previous training iteration \(k-1\) as shown in Figure 2:

\[\mathbf{G}_{L}^{k}=(1-\gamma)\frac{\partial\mathcal{L}}{\partial \mathbf{X}_{\text{in}}^{k-1}}+\gamma\frac{\partial\mathcal{L}}{\partial \mathbf{X}_{L}^{k}} \tag{13}\] \[\mathbf{G}_{l}^{k}=(1-\alpha)\tilde{\mathbf{A}}\mathbf{G}_{l+1}+ \alpha\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}^{k}}\leavevmode \nobreak\ \forall l=L-1,\dots,0 \tag{14}\]

where \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}^{k}}\) is approximated by \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}^{k}}\), and \(\mathbf{G}_{0}\) provides an approximation for gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}^{k}}\) so that the gradient of model parameters, i.e., \(\frac{\partial\mathcal{L}}{\partial\mathbf{G}^{k}}\), can be further computed by chain rules. Similar to the forward computation, the momentum correction in Eq. (13) compensates the gradient changes over training iterations by mixing the current gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{L}^{k}}\) and history gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}^{k-1}}\) with hyperparameter \(\gamma\).

**Remark 3** (Computation and memory efficiency).: _Since LazyGNN uses very few aggregation layers, it significantly reduces the computation and memory cost. Moreover, LazyGNN does not need to store intermediate hidden activation values in the aggregation layers because the computation of implicit gradient is agnostic to the forward propagation trajectory, which further reduces memory cost._

### Mini-batch LazyGNN with Subgraph Sampling

As illustrated in Section 3.1, LazyGNN solves the inherent neighborhood explosion problem so that when handling large-scale graphs, the mini-batch sampling only needs to sample neighboring nodes within a few hop distances. To further demonstrate the superior scalability, we introduce a mini-batch LazyGNN to further improve the computation and memory efficiency with efficient data sampling as shown in Figure 3. In each training iteration \(k\), we sample a target node set \(V_{1}^{k}\) and their \(L\)-hop neighbor set \(V_{2}^{k}\), and we denote the union of these two nodes as \(V^{k}=V_{1}^{k}\cup V_{2}^{k}\). An induced subgraph \(\tilde{\mathbf{A}}_{V^{k}}\) is constructed based on the node set \(V^{k}\). Note that LazyGNN works well with small \(L\in\{1,2\}\) so that the data sampling is extremely efficient.

**Forward Computation.** The forward computation of mini-batch LazyGNN works as follows:
\[(\mathbf{X}_{\text{in}}^{k})_{V_{k}} =f\big{(}(\mathbf{X}_{\text{fea}})_{V^{k}},\Theta^{k}\big{)}, \tag{15}\] \[(\mathbf{X}_{0})_{V_{k}} =(1-\beta)(\mathbf{M}_{\text{fea}})_{V^{k}}+\beta(\mathbf{X}_{ \text{in}}^{k})_{V_{k}}\] (16) \[(\mathbf{X}_{l+1}^{k})_{V^{k}} =(1-\alpha)\tilde{\mathbf{A}}_{V^{k}}(\mathbf{X}_{l}^{k})_{V^{k}} +\alpha(\mathbf{X}_{\text{in}}^{k})_{V^{k}},\forall l\] (17) \[(\mathbf{M}_{\text{fea}})_{V_{1}^{k}} =(\mathbf{X}_{L}^{k})_{V_{1}^{k}} \tag{18}\]

The node features \((\mathbf{X}_{\text{fea}})_{V^{k}}\) for the node set \(V^{k}\) are sampled to form a mini-batch. The corresponding diffused node features \((\mathbf{M}_{\text{fea}})_{V^{k}}\) of the same node set are retrieved from the feature storage \(\mathbf{M}_{\text{fea}}\) and then combined with current node features \((\mathbf{X}_{\text{in}}^{k})_{V_{k}}\) in Eq. (16). The lazy forward propagation runs on the subgraph \(\tilde{\mathbf{A}}_{V^{k}}\) in Eq. (17). Finally, \((\mathbf{X}_{L}^{k})_{V_{1}^{k}}\) provides the prediction for target nodes \(V_{1}^{k}\), which is further maintained in the features storage \(\mathbf{M}_{\text{fea}}\).

**Backward Computation.** The backward computation of mini-batch LazyGNN works as follows:

\[(\mathbf{G}_{L}^{k})_{V^{k}}=(1-\gamma)(\mathbf{M}_{\text{grad}} )_{V^{k}}+\gamma\frac{\partial\mathcal{L}}{\partial(\mathbf{X}_{L}^{k})_{V^{k}}} \tag{19}\] \[(\mathbf{G}_{l}^{k})_{V^{k}}=(1-\alpha)\tilde{\mathbf{A}}_{V^{k} }(\mathbf{G}_{l+1})_{V^{k}}+\frac{\alpha}{\partial(\mathbf{X}_{L}^{k})_{V^{k}} }\;\forall l\] (20) \[(\mathbf{M}_{\text{grad}})_{V_{1}^{k}}=(\mathbf{G}_{0}^{k})_{V_{1 }^{k}} \tag{21}\]

In Eq. (19), \((\mathbf{M}_{\text{grad}})_{V^{k}}\), the history gradient with respect to node features \((\mathbf{X}_{\text{in}}^{k-1})_{V^{k}}\), is retrieved from the gradient storage \(\mathbf{M}_{\text{grad}}\) and then combined with \(\frac{\partial\mathcal{L}}{\partial(\mathbf{X}_{L}^{k})_{V^{k}}}\), the gradient with respect to current node features \((\mathbf{X}_{L}^{k})_{V^{k}}\). The lazy backward propagation runs on the subgraph \(\mathbf{A}_{V^{k}}\) in Eq. (20). Finally, \((\mathbf{G}_{0}^{k})_{V_{1}^{k}}\) provides an estimation for \(\frac{\partial\mathcal{L}}{\partial(\mathbf{X}_{\text{in}}^{k})_{V_{1}^{k}}}\), which is used to compute gradient of parameters \(\Theta^{k}\) by chain rules and then maintained in the gradient storage \(\mathbf{M}_{\text{grad}}\).

### Complexity Analysis

As discussed in Section 3.1 and Section 3.2, LazyGNN is scalable due to its efficiency in computation, memory, and data sampling. In this section, we provide complexity analyses to further demonstrate its superior scalability. Since the major efficiency difference lies in feature aggregation layers, we focus on the complexity of feature aggregation. Suppose \(L\) is the number of propagation layers, \(H\) is the size of hidden units, \(N\) is the number of nodes, \(M\) is the number of edges. For simplicity, we assume that \(H\) is the same for all hidden layers. Performing one feature aggregation requires a sparse-dense matrix multiplication that needs \(\mathcal{O}(MH)\) operations. Therefore, the computation complexity for forward feature aggregations and backward gradient aggregations is \(\mathcal{O}(2LMH)\) per epoch. However, the number of layers \(L\) in LazyGNN is much smaller than existing approaches, so it significantly reduces the computation cost.

For memory complexity, \(\mathcal{O}(NH)\) space is required for the storage of history feature \(\mathbf{X}_{L}^{k}\) and gradient \(\frac{\partial\mathcal{L}}{\partial\mathbf{X}_{\text{in}}^{k}}\) in LazyGNN. LazyGNN does not store the intermediate state at each feature aggregation layer because the backward gradient computation is agnostic to the forward computation trajectory. Therefore, the total memory complexity for LazyGNN is \(\mathcal{O}(NH)\), which is independent of the number of aggregation layers. In contrast, existing state-of-art history-embedding based models, such as GNNAutoScale (Fey et al., 2021) and GraphFM (Yu et al., 2022), require the storage of feature embeddings in every layer, which leads to the memory complexity \(\mathcal{O}(LNH)\). Moreover, they also require the storage of all intermediate layers for gradient computation. Their memory cost increases linearly with the number of layers, which is essential in capturing long-distance dependency in those models. In fact, because of the large storage requirements, GAS chooses to save the history embedding in CPU memory, but the data movement between CPU and GPU can be dominant compared with the movement in GPU as verified in Section 4.2.

Mini-batch LazyGNN further reduces the computation and memory cost by data sampling, which enjoys the same benefits as existing sampling-based methods. Compared to classic sampling-based models, such as GraphSAGE and GraphSAINT, mini-batch LazyGNN only needs to sample neighbors in one hop or two hops distance so that the data sampling is efficient. Above analyses indicate that LazyGNN is highly efficient in computation, memory, and data sampling. The practical running time and memory cost will be discussed in Section 4.2. Note that although we do not investigate the distributed setting in this paper, LazyGNN is promising for its parallelism since it requires fewer communication rounds due to lazy propagation.

## 4 Experiments

In this section, we provide comprehensive experiments to validate the superior prediction performance and scalability of LazyGNN. Specifically, we try to answer the following questions: (Q1) Can LazyGNN achieve strong prediction accuracy on large-scale graph benchmarks? (Section 4.1) (Q2) Can LazyGNN handle large graphs more efficiently than existing approaches? (Section 4.2) (Q3) What are the impacts of the hyperparameters in LazyGNN? (Section 4.3)

Figure 3: Mini-batch LazyGNN with Feature & Gradient Storage.
### Prediction Performance

**Experimental settings.** We conduct experiments on multiple large-scale graph datasets including REDDIT, YELP, FLICKR, ogbn-arxiv, and ogbn-products (Hu et al., 2020). We evaluate the graph representation learning by node classification accuracy in the semi-supervised setting. We provide a performance comparison with multiple baselines including GCN (Kipf and Welling, 2016), GraphSAGE (Hamilton et al., 2017), FastGCN (Chen et al., 2018), LADIES (Zou et al., 2019), VR-GCN (Chen et al., 2017), MVS-GNN (Cong et al., 2020), Cluster- GCN (Chiang et al., 2019), GraphSAINT (Zeng et al., 2020), SGC (Wu et al., 2019), SIGN (Rossi et al., 2020), GAS (Fey et al., 2021). The hyperparameter tuning of baselines closely follows the setting in GNAAutoScale (Fey et al., 2021). For LazyGNN, we set \(\beta=0.5\), and \(\gamma=0.5\). Other hyperparameters are tuned from the following search space: (1) learning rate: \(\{0.01,0.001,0.0001\}\); (2) weight decay: \(\{0.5e-4,5e-5\}\); (3) dropout: \(\{0.1,0.3,0.5,0.7\}\); (4) propagation layers : \(L\in\{1,2\}\); (5) MLP layers: \(\{3,4\}\); (6) MLP hidden units: \(\{256,512\}\); (7) \(\alpha\in\{0.01,0.1,0.2,0.5,0.8\}\). For data sampling, while LazyGNN is compatible with any sampling method, we adopt the subgraph sampling strategy in GNAutoScale to ensure a fair comparison.

**Performance analysis.** From the prediction accuracy as summarized in Table 1, we can make the following observations: (1) LazyGNN achieves state-of-art performance on almost all datasets. In particular, LazyGNN achieves \(73.2\%\) and \(81.0\%\) accuracy on ogbn-arxiv and ogbn-products. The only exceptions are that GraphSAINT slightly outperforms LazyGNN by \(0.2\%\) on REDDIT and that GCNII outperforms LazyGNN by \(1.1\%\) on FLICKR. However, LazyGNN is much more efficient than GraphSAINT and GCNII in computation and memory cost. Note that GCNII achieves the reported performance with 10 transformation layers and 8 propagation layers while LazyGNN only uses 4-layer MLP and 2-layer propagation. We believe a stronger performance can be achieved with a more thorough architecture tuning on LazyGNN; (2) Importantly, we observe that full-batch LazyGNN closely matches mini-batch LazyGNN, which indicates that the subgraph sampling in LazyGNN can significantly improve the scalibility (as will be discussed in Section 4.2) but do not heavily hurt the prediction performance. Sometimes, LazyGNN even improves the generalization performance due to sampling (such as for REDDIT); (3) LazyGNN also consistently outperforms APPNP with 10 propagation layers, which indicates the advantage of capturing long-distance dependency in graphs by lazy propagation. The reproduction of GAS+GCN on YELP is much worse than reported, therefore we omit the result.

**OGB leaderboard.** To further demonstrate the superior prediction performance of LazyGNN, we also compare LazyGNN with the best-ranked GNN models on the OGB Leaderboard as of Jan. 26, 20231. To make a fair comparison, we only consider the models in the leaderboard that _do not use any external data_ such as the raw text feature. For ogbn-arxiv dataset, we choose AGDN (Sun and Wu, 2020), LGGNN (Ma, 2022), C&S (Huang et al., 2020), EnGCN (Duan et al., 2022) as baselines. For ogbn-products dataset, we compare with SCR (Zhang et al., 2021), SAGN (Sun et al., 2021), GAMLP (Zhang et al., 2022), C&S (Huang et al., 2020), and EnGCN (Duan et al., 2022). Note that these top-ranked models use additional techniques to boost the performance. For instance, AGDN uses Bags of Tricks (BoT) (Wang et al., 2021), self-Knowledge Distillation (self-KD) (Sun and Wu, 2020), and Correct and Smooth (C&S) (Huang et al., 2020); LGGNN uses C&S and Label Reuse (Ma, 2022); GAMLP uses Reliable Label Utilization (RLU) (Zhang et al., 2022) and SCR (Zhang et al., 2021); EnGCN uses self-label enhancement (SLE) (Sun et al., 2021), Label Propagation (LP), and Ensemble with weighted majority voting (MV). LazyGNN is also compatible with these techniques, and we only integrate self label enhancement (SLE) (Sun et al., 2021) to improve LazyGNN.

Footnote 1: [https://ogb.stanford.edu/docs/leader_nodeprop/](https://ogb.stanford.edu/docs/leader_nodeprop/)

The performance on ogbn-arxiv and ogbn-products is summarized in Table 2 and Table 3, respectively. They clearly demonstrate that LazyGNN (mini-batch) achieves the state-of-art performance compared with these top-ranked GNN models on the highly competitive leaderboard. LazyGNN significantly outperforms many strong baselines such as

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline
**\# nodes** & 230K & 89K & 717K & 169K & 2.4M \\
**\# edges** & 11.6M & 450K & 7.9M & 1.2M & 61.9M \\
**Method** & Reddit & Flickr & Yelp & ogbn- & ogbn- \\  & & & & & & \\ \hline GraphSAGE & 95.4 & 50.1 & 63.4 & 71.5 & 78.7 \\ FastGCN & 93.7 & 50.4 & — & — & — \\ LADIES & 92.8 & — & — & — & — \\ VR-GCN & 94.5 & — & 61.5 & — & — \\ MVS-GNN & 94.9 & — & 62.0 & — & — \\ Cluster-GCN & 96.6 & 48.1 & 60.9 & — & 79.0 \\ GraphSAINT & **97.0** & 51.1 & 65.3 & — & 79.1 \\ SGC & 96.4 & 48.2 & 64.0 & — & — \\ SIGN & 96.8 & 51.4 & 63.1 & — & 77.6 \\ \hline GCN (full) & 95.4 & 53.7 & OOM & 71.6 & OOM \\ APPNP (full) & 96.1 & 53.4 & OOM & 71.8 & OOM \\ GCNII (full) & 96.1 & 55.3 & OOM & 72.8 & OOM \\ \hline GCN (GAS) & 95.4 & 54.0 & — & 71.7 & 76.7 \\ APPNP (GAS) & 96.0 & 52.4 & 63.8 & 71.9 & 76.2 \\ GCNII (GAS) & 96.7 & **55.3** & 65.1 & 72.5 & 77.2 \\ \hline LazyGNN (full) & 96.2 & 54.2 & OOM & **73.2** & OOM \\ LazyGNN (mini) & 96.8 & 54.0 & **65.4** & 73.0 & **82.0** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Prediction accuracy (\(\%\)) on large-scale graph datasets. “full” and “mini” stand for full-batch and mini-batch respectively. OOM stands for out-of-memory.
AGDN, LGGNN, GAMLP, and SAGN, and it slightly outperforms the best-performing EnGCN (Ensemble GCN) with significantly better efficiency (1.5 hours vs 2.4 hours).

### Efficiency Analysis

To verify the scalability of LazyGNN, we provide empirical efficient analysis compared with state-of-art end-to-end trainable and scalable GNNs such as GNAAutoScale (GAS) since it has been shown to be the most efficient algorithm for training large-scale GNNs (Fey et al., 2021). Specifically, we measure the memory usage and running time for the training on ogbn-products dataset using the same batch size. We run GAS using the authors' code which stores the feature embedding for each layer in CPU memory. For LazyGNN, the memory for storing two small tensors (i.e., feature and gradient) is required. Although the storage memory for LazyGNN can easily fit into the GPU, we measure the running time for two cases depending on whether the storage is on CPU or GPU memory. Note that all methods use similar number of training epochs as will be verified in Section 4.3, so we compare them using the running time per epoch. For GAS baselines, we use the architecture that can reach the best performance.

From the measurements summarized in Table 4, we can make the following observations: (1) LazyGNN (GPU) has a much shorter running time (7.7s) per epoch compared with all baselines; GAS becomes much slower due to the memory movement between GPU and CPU. (2) LazyGNN (GPU+CPU) is faster than APPNP (GAS) because APPNP requires more data movement for feature propagation. This further verifies the computation and memory efficiency of LazyGNN. (3) LazyGNN is memory efficient since only 878 MB additional storage is required to store two small tensors for feature and gradient regardless of the number of layers. In contrast, GAS requires large memory storage due to the feature storage for each layer. These observations verify that LazyGNN is efficient in memory and computation.

### Ablation Study

We provide detailed ablation studies on the impacts of hyperparameter settings in LazyGNN using ogbn-arixv dataset.

**Lazy Propagation.** We conduct the ablation study to analyze the impact of feature propagation layers. Specifically, we use the same MLP for LazyGNN and APPNP, and we change the number of propagation layers \(K\). The comparison in Table 5 demonstrates that: (1) LazyGNN have much better performance than APPNP when using the same number of propagation layers; (2) LazyGNN archives great performance with even \(1\) propagation layer (\(72.5\%\)), and it already arrives at the best performance (\(73.0\%\)) with \(2\) propagation layers. These results confirm LazyGNN's advantages and capability in capturing long-distance dependency in graphs with very few feature propagation layers.

**Sensitivity Analysis.** We provide a detailed sensitivity analysis for the two additional hyperparameters \(\beta\) and \(\gamma\) in LazyGNN. The results in Table 6 show that: (1) \(\{\beta=0,\gamma=0\}\) produces the worst performance (\(71.0\%\)) because it fully trusts the history embedding that might be outdated due to the additional variations caused by dropout as demonstrated in Section 2.2; (2) \(\{\beta=1,\gamma=1\}\) performs slightly better (\(71.7\%\)) even without lazy propagation because 2 propagation layers can provide an approximate solution; (3) \(\{\beta=0.5,\gamma=0.5\}\) performs the best (\(73.0\%\)) because it achieves a good trade-off between historical information and current iteration. It compensates for the staleness in the history storage as designed; (4) The performance of LazyGNN is quite stable for a large range setting of \(\beta\) and \(\gamma\)

\begin{table}
\begin{tabular}{l c c} \hline \hline
**Method** & \begin{tabular}{c} GPU \\ Memory \\ \end{tabular} & \begin{tabular}{c} CPU \\ Memory \\ \end{tabular} & 
\begin{tabular}{c} Running \\ Time \\ \end{tabular} \\ \hline GCN (GAS) & 2512 & 4783 & 28.8 \\ GCNII (GAS) & 3035 & 4783 & 44.3 \\ APPNP (GAS) & 2142 & 3952 & 84.5 \\ LazyGNN (GPU+CPU) & 2101 & 878 & 44.5 \\ LazyGNN (GPU) & 2981 & — & 7.7 \\ \hline \hline \end{tabular}
\end{table}
Table 4: Memory usage (MB) and running time for each epoch (seconds per epoch) on ogbn-products.

\begin{table}
\begin{tabular}{l c} \hline \hline
**Method** & 
\begin{tabular}{c} Accuracy (\%) \\ \end{tabular} \\ \hline GAMLP + RLU + SCR & 85.05 \\ GAMLP + RLU + SCR + C\&S & 85.20 \\ SAGN + SLE & 84.68 \\ SAGN + SLE + C\&S & 84.85 \\ EnGCN + SLE + LP + MV & 87.98 \\ LazyGNN + SLE & **88.05** \\ \hline \hline \end{tabular}
\end{table}
Table 3: Comparison with top-ranked models on ogbn-products.

\begin{table}
\begin{tabular}{l c} \hline \hline
**Method** & 
\begin{tabular}{c} Accuracy (\%) \\ \end{tabular} \\ \hline APPNP (K=1) & 70.6 \\ APPNP (K=2) & 71.4 \\ APPNP (K=3) & 71.5 \\ LazyGNN (K=1) & 72.5 \\ LazyGNN (K=2) & 73.0 \\ LazyGNN (K=3) & 73.0 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Prediction accuracy with different numbers of feature propagation layers on ogbn-arxiv.
In fact, we fix \(\beta=\gamma=0.5\) for LazyGNN in all experiments for simplicity. Therefore, it requires almost negligible effort for additional hyperparameter tuning.

**Convergence.** We show the convergence of LazyGNN during the training process using ogbn-arxiv dataset. The convergence of validation accuracy in Figure 4 demonstrates that LazyGNN has a comparable convergence speed with GCN (GAS) and GCNII (GAS), and is slightly faster than APPNP (GAS) in terms of the number of training epochs.

## 5 Related Work

It has been generally demonstrated that it is beneficial to capture long-distance relations in graphs by stacking more feature aggregation layers or unrolling various fixed point iterations in GNNs (Gasteiger et al., 2018; Gu et al., 2020; Liu et al., 2020; Chen et al., 2020; Li et al., 2021; Ma et al., 2020; Pan et al., 2020; Zhu et al., 2021; Chen et al., 2020). But these works suffer from scalability concerns due to the neighborhood explosion problem (Hamilton et al., 2017).

A large body of existing research focuses on improving the efficiency and scalability of large-scale GNNs using various novel designs, such as sampling methods, pre-computing or post-computing methods, and distributed methods. Sampling methods adopt mini-batch training strategies to reduce computation and memory requirements by sampling nodes and edges. They mitigate the neighbor explosion issue by either removing neighbors (Hamilton et al., 2017; Chen et al., 2018; Zeng et al., 2020; Zou et al., 2019) or updating with feature memory (Fey et al., 2021; Yu et al., 2022). Pre-computing or post-computing methods separate the feature aggregation and prediction models into two stages, such as pre-computing the feature aggregation before training (Wu et al., 2019; Rossi et al., 2020; Sun et al., 2021; Zhang et al., 2022; Bojchevski et al., 2020) or post-processing with label propagation after training (Zhu, 2005; Huang et al., 2020). Distributed methods distribute large graphs to multiple servers and parallelize GNNs training (Chiang et al., 2019; Chai et al., 2022; Shao et al., 2022). In retrospect, these existing approaches still suffer from various limitations, such as high costs in computation, memory, and communication as well as performance degradation due to large approximation errors or multi-stage training.

Different from existing works, in this work, we propose LazyGNN from a substantially different and novel perspective and propose to capture the long-distance dependency in graphs by shallower models instead of deeper models. This leads to a much more efficient LazyGNN for graph representation learning. Moreover, existing approaches for scalable GNNs can be used to further accelerate LazyGNN. The proposed mini-batch LazyGNN is a promising example.

LazyGNN also draws insight from existing research on graph signal processing and implicit modeling. The forward construction of LazyGNN follows the optimization perspective of GNNs (Ma et al., 2020), and it is can be easily generalized to other fixed point iterations with various diffusion properties. LazyGNN is related to recent works in implicit modeling, such as Neural Ordinary Differential Equations (Chen et al., 2018), Implicit Deep Learning (El Ghaoui et al., 2021), Deep Equilibrium Models (Bai et al., 2019), and Implicit GNNs (Gu et al., 2020). LazyGNN focuses on developing shallow models with highly scalable and efficient computation while these implicit models demand significantly more computation resources.

## 6 Conclusions

In this paper, we propose LazyGNN, a novel shallow model to solve the neighborhood explosion problem in large-scale GNNs while capturing long-distance dependency in graphs through lazy propagation. We also develop a highly scalable and efficient variant, mini-batch LazyGNN, to handle large graphs. Comprehensive experiments demonstrate its superior prediction performance and efficiency on large-scale graph datasets. We plan to explore lazy propagation for other types of GNN models and its further acceleration based on other scalable techniques in the future.

\begin{table}
\begin{tabular}{c c} \hline \hline
**Hyperparameter settings** & Accuracy (\%) \\ \hline \(\beta=0.0,\;\gamma=0.0\) & 71.0 \\ \(\beta=1.0,\;\gamma=1.0\) & 71.7 \\ \hline \(\beta=0.5,\;\gamma=0.1\) & 72.5 \\ \(\beta=0.5,\;\gamma=0.5\) & 73.0 \\ \(\beta=0.5,\;\gamma=0.9\) & 72.8 \\ \hline \(\beta=0.1,\;\gamma=0.5\) & 72.4 \\ \(\beta=0.5,\;\gamma=0.5\) & 73.0 \\ \(\beta=0.9,\;\gamma=0.5\) & 72.7 \\ \hline \hline \end{tabular}
\end{table}
Table 6: Prediction accuracy of mini-batch LazyGNN for different hyperparameter settings on ogbn-arxiv.

Figure 4: Validation performance versus training epochs.