# On the Generalization Capability of Temporal Graph Learning Algorithms:
Theoretical Insights and a Simpler Method

Weilin Cong  Penn State University  Jian Kang  University of Rochester
Hanghang Tong  University of Illinois at Urbana-Champaign  Mehrdad Mahdavi
Penn State University

###### Abstract

Temporal Graph Learning (TGL) has become a prevalent technique across diverse
real-world applications, especially in domains where data can be represented
as a graph and evolves over time. Although TGL has recently seen notable
progress in algorithmic solutions, its theoretical foundations remain largely
unexplored. This paper aims at bridging this gap by investigating the
generalization ability of different TGL algorithms (e.g., GNN-based, RNN-
based, and memory-based methods) under the finite-wide over-parameterized
regime. We establish the connection between the generalization error of TGL
algorithms and 1 “the number of layers/steps” in the GNN-/RNN-based TGL
methods and 2 “the feature-label alignment (FLA) score”, where FLA can be used
as a proxy for the expressive power and explains the performance of memory-
based methods. Guided by our theoretical analysis, we propose Simplified-
Temporal-Graph-Network (SToNe), which enjoys a small generalization error,
improved overall performance, and lower model complexity. Extensive
experiments on real-world datasets demonstrate the effectiveness of SToNe. Our
theoretical findings and proposed algorithm offer essential insights into TGL
from a theoretical standpoint, laying the groundwork for the designing
practical TGL algorithms in future studies.

##  1 Introduction

Temporal graph learning (TGL) has emerged as an important machine learning
problem and is widely used in a number of real-world applications, such as
traffic prediction [Yuan and Li, 2021, Zhang et al., 2021], knowledge graphs
[Cai et al., 2022, Leblay and Chekol, 2018], and recommender systems [Kumar et
al., 2019, Rossi et al., 2020, Xu et al., 2020a]. A typical downstream task of
temporal graph learning is link prediction, which focuses on predicting future
interactions among nodes. For example in an online video recommender system,
the user-video clicks can be modeled as a temporal graph whose nodes represent
users and videos, and links are associated with timestamps indicating when
users click videos. Link prediction between nodes can be used to predict if
and when a user is interested in a video. Therefore, designing graph learning
models that can capture node evolutionary patterns and accurately predict
future links is important.

TGL is generally more challenging than static graph learning, thereby
requiring more sophisticated algorithms to model the temporal evolutionary
patterns [Huang et al., 2023]. In recent years, many TGL algorithms [Kumar et
al., 2019, Xu et al., 2020a, Rossi et al., 2020, Sankar et al., 2020, Wang et
al., 2021e] have been proposed that leverage memory blocks, self-attention,
time-encoding function, recurrent neural networks (RNNs), temporal walks, and
message passing to better capture the meaningful structural or temporal
patterns. For instance, JODIE [Kumar et al., 2019] maintains a memory block
for each node and utilizes an RNN to update the memory blocks upon the
occurance of each interaction; TGAT [Xu et al., 2020a] utilizes self-attention
message passing to aggregate neighbor information on the temporal graph; TGN
[Rossi et al., 2020] combines memory blocks with message passing to allow each
node in the temporal graph to have a receptive field that is not limited by
the number of message-passing layers; DySAT [Sankar et al., 2020] uses self-
attention to capture structural information and uses RNN to capture temporal
dependencies; CAW [Wang et al., 2021e] captures temporal dependencies between
nodes by performing multiple temporal walks from the root node and
representing neighbors’ identities via the probability of a node appearing in
the temporal walks. GraphMixer [Cong et al., 2023] combines MLP-mixer
[Tolstikhin et al., 2021] with 1-hop most recent neighbor aggregation for
temporal graph learning.

Figure 1: Relationship between the generalization error (in Theorem 1) and
empirical average precision score. Generalization error (GE) and average
precision (AP) have an inverse correlation, i.e., the larger the GE, the lower
the AP. Each marker is one experiment run. The same method’s GE changes at
each run because it depends on feature-label alignment, which changes with
different weight initialization. More details on the computation of GE are
deferred to Appendix I.4.

Although empirically powerful, the theoretical foundations of the
aforementioned TGL algorithms remain largely under-explored. A recent study
[Cong et al., 2023] have demonstrated that simple models could still exhibit
superior performance compared to more complex models [Kumar et al., 2019, Xu
et al., 2020a, Rossi et al., 2020, Sankar et al., 2020, Wang et al., 2021e],
but rigorous understanding of this phenomena remains elusive. Therefore, it is
imperative to unravel the learning capabilities of these methods from the
theoretical perspective, which will be in turn beneficial to design a more
practically efficacious TGL algorithm in the future studies. Recently, [Souza
et al., 2022, Gao and Ribeiro, 2022] study the expressive power of temporal
graphs by extending the 1-WL (Weisfeiler-Lehman) test to the temporal graph.
[Souza et al., 2022] shows that using an injective aggregation function in TGL
models is necessary to be most expressive, and that increasing the number of
layers/steps in GNN-/RNN-based TGL methods or using memory blocks can increase
the WL-based expressive power. However, WL-based expressive power analysis may
not 1 explain why one TGL method performs better than another when neither of
them uses injective aggregation functions (e.g., mean-average aggregation), 2
explain why existing TGL methods prefer to use shallow GNN structure (e.g.,
[Xu et al., 2020a, Fan et al., 2021, Rossi et al., 2020] uses 1- or 2-layer in
official implementation) or a small number of RNN steps (e.g., [Wang et al.,
2021e, Sankar et al., 2020] use less than 4 steps in official implementation)
despite the theory suggesting deeper/larger structures, and 3 explain how well
TGL algorithms generalize to the unseen data.

As an alternative, we directly study the generalization ability of different
TGL algorithms. We establish the generalization error bound of GNN-based [Xu
et al., 2020a, Fan et al., 2021], RNN-based [Wang et al., 2021e], and memory-
based [Kumar et al., 2019] TGL algorithms under the finite-wide over-
parameterized regime in Theorem 1 via deep learning theory. In particular, we
show that the generalization error bound decreases with respect to the number
of training data, but increases with respect to the number of layers/steps in
the GNN-/RNN-based methods and the feature-label alignment (FLA), where FLA is
defined as the projection of labels on the inverse of empirical neural tangent
kernel
$\mathbf{y}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y},\leavevmode\nobreak\
\mathbf{J}=[\text{vec}(\nabla_{\theta}f_{i}(\bm{\theta}_{0}))]_{i=1}^{N}$ as
in Definition 1. We show that the FLA could potentially 1 serve as a proxy for
the expressive power measurement, 2 explain the impact of input data selection
on model performance, and 3 explain why the memory-based methods do not
outperform GNN-/RNN-based methods, even though their generalization error is
less affected as the number of layers/steps increases. Our generalization
bound demonstrates the importance of “selecting proper input data” and “using
simpler model architecture”, which have been previously empirically observed
in [Cong et al., 2023], but lack theoretical understanding.

Guided by our theoretical analysis, we propose Simplified-Temporal-Graph-
Network (SToNe), which achieves compatible performance on most real-world
datasets, but with less weight parameters and lower computational cost. As
shown in Figure 1, our proposed SToNe enjoys not only small generalization
error theoretically but also a compatible overall average precision score.
Extensive ablation studies are conducted on real-world datasets to demonstrate
the effectiveness of SToNe.

#### Contributions.

We summarize the main contributions as follows:

  * •

(_Theory_) In Section 3, we analyze the generalization error of memory-based,
GNN-based, and RNN-based methods in Theorem 1. Then, in Section 3.3, we reveal
the relationship between generalization error and the number of layers/steps
in the GNN-/RNN-based method and feature-label alignment (FLA), where FLA
could also be potentially used as a proxy for the expressive power and explain
the performance of the memory-based method.

  * •

(_Algorithm_) In Section 4, guided by our theoretical analysis, we propose
SToNe which not only enjoys a small theoretical generalization error but also
bears a better overall empirical performance, a smaller model complexity, and
a simpler architecture compared to baseline methods.

  * •

(_Experiments_) In Section 5 and Appendix B, we conduct extensive experiments
and ablation studies on real-world datasets to demonstrate the effectiveness
of our proposal SToNe .

##  2 Related works and preliminaries

In this section, we briefly summarize related works and preliminaries. The
more detailed discussions on expressive power, generalization, and existing
TGL algorithms are deferred to Appendix H.

Generalization and expressive power. Statistical learning theories have been
used to study the generalization of GNNs, including uniform stability [Verma
and Zhang, 2019, Zhou and Wang, 2021, Cong et al., 2021], Rademacher
complexity [Garg et al., 2020, Oono and Suzuki, 2020, Du et al., 2019], PAC-
Bayesian [Liao et al., 2020], PAC-learning [Xu et al., 2020c], and uniform
convergence [Maskey et al., 2022]. Besides, [Souza et al., 2022, Gao and
Ribeiro, 2022] study the expressive power of temporal graph networks via graph
isomorphism test. However, the aforementioned analyses are data-independent
and only dependent on the number of layers or hidden dimensions, which cannot
fully explain the performance difference between different algorithms. As an
alternative, we study the generalization error of different TGL methods under
the finite-wide over-parameterized regime [Xu et al., 2021a, Arora et al.,
2019, Arora et al., 2022], which not only is model-dependent but also could
capture the data dependency via feature-label alignment (FLA) in Definition 1.
Most importantly, FLA could be empirically computed, reveals the impact of
input data selection on model performance, and could be used as a proxy for
the expressive power of different algorithms.

Figure 2: An illustration of temporal graph data with nodes
$v_{1},\ldots,v_{5}$ and timestamps $t_{1},\ldots,t_{6}$ that indicate when
two nodes interact.

Temporal graph learning. Figure 2 is an illustration of a temporal graph,
where each node has node feature $\mathbf{x}_{i}$, each node pair could have
multiple temporal edges with different timestamps $t$ and edge features
$\mathbf{e}_{ij}(t)$. We classify several chosen representative temporal graph
learning methods into memory-based (e.g., JODIE), GNN-based (e.g., TGAT,
TGSRec), memory&GNN-based (e.g., TGN, APAN, and PINT), RNN-based (e.g., CAW),
and GNN&RNN-based (e.g., DySAT) methods. For example, JODIE [Kumar et al.,
2019] maintains a memory block for each node and updates the memory block by
an RNN upon the happening of each interaction; TGAT [Xu et al., 2020a] first
constructs the temporal computation graph, then recursively computes the
hidden representation of each node using GNN; TGN [Rossi et al., 2020] first
uses memory blocks to capture all temporal interactions via JODIE then applies
TGAT on the latest representation of the memory blocks of each node to capture
the spatial information; CAW [Wang et al., 2021e] proposes to first construct
a set of sequential temporal events via random walks, then use RNN to
aggregate the information from temporal events; DySAT [Sankar et al., 2020]
uses GNN to extract the spatial features and applies RNN on the output of GNN
to capture the temporal dependencies.

##  3 Generalization of temporal graph learning methods

Firstly, we introduce the problem setting of theoretical analysis in Section
3.1. Then, we formally define the feature-label alignment (FLA) and derive the
generalization bound in Section 3.2. Finally, we provide discussion on the
generalization bound and its connection to FLA in Section 3.3.

###  3.1 Problem setting for theoretical analysis

For theoretical analysis, let us suppose the TGL models receives a sequence of
data points $(X_{1},y_{1}),...,(X_{N},y_{N})$ that is the realization of some
non-stationary process, where $X_{i}$ and $y_{i}$ stand for a single input
data (e.g., a root node and its induced subgraph or random-walk path) and its
label at the $i$-th iteration. During training, stochastic gradient descent
(SGD) first uses $(X_{1},y_{1})$ and the initial model $f(\bm{\theta}_{0})$ to
generate the $\bm{\theta}_{1}$. Next, SGD uses the second example
$(X_{2},y_{2})$ and the previously obtained model $f(\bm{\theta}_{1})$ to
generate $\bm{\theta}_{2}$, and so on. The training process is outlined in
Algorithm 1. We aim to develop a unified theoretical framework to examine the
generalization capability of the three most fundamental types of TGL methods,
including GNN-based method, RNN-based method, and memory-based method on a
node-level binary classification task. Our goal is to upper bound the expected
0-1 error
$\mathbb{E}[\text{loss}^{0-1}_{N}(\widetilde{\bm{\theta}})|\mathcal{D}_{1}^{N-1}]$
conditioned on a sequence of sampled data points $\mathcal{D}_{1}^{N-1}$,
where
$\text{loss}_{i}^{0-1}(\bm{\theta})=\mathbf{1}\\{y_{i}f_{i}(\bm{\theta})<0\\}$
is the 0-1 loss computed by model $f(\bm{\theta})$ on data $(X_{i},y_{i})$ at
$i$-th iteration, and
$\mathcal{D}_{1}^{i-1}=\\{(X_{1},y_{1}),\ldots,(X_{i-1},y_{i-1})\\}$ is the
set of all data points sampled before the $i$-th iteration. Our analysis of
the fundamental TGL methods can pave the way to understand more advanced
algorithms, such as GNN&RNN-based and memory&GNN-based methods.

Algorithm 1 Using SGD for temporal graph learning

Initialize $\bm{\theta}_{0}$ from Gaussian distribution $\mathcal{N}(0,1/m)$
with $m$ as hidden dimension.

for $i=1,2,\ldots,N$ do

Sample an input data $(X_{i},y_{i})$ condition on all the previous sampled
data $\mathcal{D}_{1}^{i-1}$

Predict the label of $X_{i}$ as $f_{i}(\bm{\theta}_{i-1})$

Update weight parameters by
$\bm{\theta}_{i}=\bm{\theta}_{i-1}-\eta\nabla_{\theta}\leavevmode\nobreak\
\text{loss}(f_{i}(\bm{\theta}_{i-1}),y_{i})$.

end for

Return weight parameters by selecting $\widetilde{\bm{\theta}}$ randomly from
$\\{\bm{\theta}_{0},\ldots,\bm{\theta}_{N-1}\\}$.

GNN-based method. We compute the representation of node $v_{i}$ at time $t$ by
applying GNN on the temporal graph that only considers the temporal edges with
timestamp before time $t$. The GNN has
$\bm{\theta}=\\{\mathbf{W}^{(\ell)}\\}_{\ell=1}^{L}$ as the parameters to
optimize and $\alpha\in\\{0,1\\}$ as a binary variable that controls whether a
residual connection is used. The final prediction on node $v_{i}$ is computed
as $f_{i}(\bm{\theta})=\mathbf{W}^{(L)}\mathbf{h}_{i}^{(L-1)}$, where the
hidden representation $\mathbf{h}_{i}^{(L-1)}$ is computed by

$\mathbf{h}_{i}^{(\ell)}=\sigma\big{(}\mathbf{W}^{(\ell)}\sum\nolimits_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(\ell-1)}\big{)}+\alpha\mathbf{h}_{i}^{(\ell-1)},\leavevmode\nobreak\
\mathbf{h}_{i}^{(1)}=\sigma\big{(}\mathbf{W}^{(1)}\sum\nolimits_{j\in\mathcal{N}(i)}P_{ij}\mathbf{x}_{j}\big{)}.$

Here $\sigma(\cdot)$ is the activation function, $P_{ij}$ is the aggregation
weight used for propagating information from node $v_{j}$ to node $v_{i}$, and
$\mathcal{N}(i)$ is the set of all neighbors of node $v_{i}$ in the temporal
graph. For simplicity, we assume the neighbor aggregation weights are fixed
(e.g., $P_{ij}=1/|\mathcal{N}(i)|$ for row normalized propagation) and
consider the time-encoding vector as part of the node feature vector
$\mathbf{x}_{i}$. For parameter dimension, we have
$\mathbf{W}^{(1)}\in\mathbb{R}^{m\times d}$,
$\mathbf{W}^{(\ell)}\in\mathbb{R}^{m\times m}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ 2\leq\ell\leq L-1$, and
$\mathbf{W}^{(L)}\in\mathbb{R}^{1\times m}$, where $m$ is the hidden dimension
and $d$ is the input feature dimension. TGAT [Xu et al., 2020a] can be thought
of as GNN-based method but uses self-attention neighbor aggregation.

RNN-based method. We compute the representation of node $v_{i}$ at time $t$ by
applying a multi-step RNN onto a sequence of temporal events
$\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{L-1}\\}$ that are constructed at the
target node, where each temporal event feature $\mathbf{v}_{\ell}$ is pre-
computed on the temporal graph. We consider the time-encoding vector as part
of event feature $\mathbf{v}_{\ell}$. The RNN has trainable parameters
$\bm{\theta}=\\{\mathbf{W}^{(1)},\mathbf{W}^{(2)},\mathbf{W}^{(3)}\\}$ and
$\alpha\in\\{0,1\\}$ is a binary variable that controls whether a residual
connection is used. Then, the final prediction on node $v_{i}$ is computed as
$f_{i}(\bm{\theta})=\mathbf{W}^{(3)}\mathbf{h}_{L-1}$, where the hidden
representation $\mathbf{h}_{L-1}$ is recursively compute by

$\mathbf{h}_{\ell}=\sigma\Big{(}\kappa\big{(}\mathbf{W}^{(1)}\mathbf{h}_{\ell-1}+\mathbf{W}^{(2)}\mathbf{x}_{\ell}\big{)}\Big{)}+\alpha\mathbf{h}_{\ell-1}\in\mathbb{R}^{m},$

$\sigma(\cdot)$ is the activation function, $\mathbf{h}_{0}=\mathbf{0}_{m}$ is
an all-zero vector, and $\mathbf{x}_{\ell}=\mathbf{W}^{(0)}\mathbf{v}_{\ell}$.
We normalize each $\mathbf{h}_{\ell}$ by $\kappa=1/\sqrt{2}$ so that
$\|\mathbf{h}_{\ell}\|_{2}^{2}$ does not grow exponentially with $L$. We have
$\mathbf{W}^{(1)},\mathbf{W}^{(2)}\in\mathbb{R}^{m\times m}$,
$\mathbf{W}^{(3)}\in\mathbb{R}^{1\times m}$ as trainable parameters, but
$\mathbf{W}^{(0)}\in\mathbb{R}^{m\times d}$ is non-trainable. CAW [Wang et
al., 2021e] is a special case of RNN-based method for edge classification
tasks, where temporal events are sampled by temporal walks from both the
source node and the destination node of an edge.

Memory-based method. We compute the representation of node $v_{i}$ at time $t$
by applying weight parameters on the memory block $\mathbf{s}_{i}(t)$. Let us
define $\bm{\theta}=\\{\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(4)}\\}$ as the
parameters to optimize. Then, the final prediction of node $v_{i}$ is computed
by $f_{i}(\bm{\theta})=\mathbf{W}^{(4)}\mathbf{s}_{i}(t)$ and
$\mathbf{s}_{i}(t)\in\mathbb{R}^{m}$ is updated whenever node $v_{i}$
interacts with other nodes by

$\mathbf{s}_{i}(t)=\sigma\Big{(}\kappa\big{(}\mathbf{W}^{(1)}\mathbf{s}_{i}^{+}(h_{i}^{t})+\mathbf{W}^{(2)}\mathbf{s}_{j}^{+}(h_{j}^{t})+\mathbf{W}^{(3)}\mathbf{e}_{ij}(t)\big{)}\Big{)},$

where $\sigma(\cdot)$ is the activation function, $h_{i}^{t}$ is the latest
timestamp that node $v_{i}$ interacts with other nodes before time $t$,
$\mathbf{s}_{i}(0)=\mathbf{W}^{(0)}\mathbf{x}_{i}$, and
$\mathbf{s}_{i}^{+}(t)=\text{StopGrad}(\mathbf{s}_{i}(t))$ is the memory block
of node $v_{i}$ at time $t$. We consider the time-encoding vector as part of
edge feature $\mathbf{e}_{ij}(t)$. We normalize $\mathbf{s}_{i}(t)$ by
$\kappa=1/\sqrt{3}$ so that $\|\mathbf{s}_{i}(t)\|_{2}^{2}$ does not grow
exponentially with time $t$. We have trainable parameters
$\mathbf{W}^{(1)},\mathbf{W}^{(2)}\in\mathbb{R}^{m\times m}$,
$\mathbf{W}^{(3)}\in\mathbb{R}^{m\times d}$,
$\mathbf{W}^{(4)}\in\mathbb{R}^{m}$ and fixed parameters
$\mathbf{W}^{(0)}\in\mathbb{R}^{m\times d}$. JODIE [Kumar et al., 2019] can be
thought of as memory-based method but using the attention mechanism for final
prediction.

###  3.2 Assumptions and main theoretical results

For the purpose of rigorous analysis, we make the following standard
assumptions on the feature norms [Cao and Gu, 2019, Du et al., 2018], which
could be satisfied via feature re-scaling.

######  Assumption 1.

All features has $\ell_{2}$-norm bounded by $1$, i.e., we assume
$\|\mathbf{x}_{i}\|_{2},\|\mathbf{e}_{ij}(t)\|_{2},\|\mathbf{v}_{\ell}\|\leq
1$.

In addition, we assume that the activation functions are Lipschitz continuous
in TGL methods. The following assumption holds for common activation functions
such as ReLU and LeakyReLU (which are often used in GNNs), Sigmoid and Tanh
(which are often used in RNNs).

######  Assumption 2.

The activation function has Lipschitz constant $\rho\geq 1$.

Furthermore, we make the following assumptions on the propagation matrix of
GNN models, which are previously used in [Liao et al., 2020, Cong et al.,
2021]. In practice, we know that $\tau=1$ holds for row normalized and
$\tau=\sqrt{\max_{i\in\mathcal{V}}d_{i}}$ holds for symmetrically normalized
propagation matrix.

######  Assumption 3.

The row-wise sum is bounded by
$\tau=\max_{i\in\mathcal{V}}(\sum_{j\in\mathcal{N}(i)}P_{ij})$ where $\tau\geq
1$.

Finally, we adopt the following assumption regarding the non-stationary data
generation process, which is standard assumption in time series prediction and
online learning analysis [Kuznetsov and Mohri, 2015, Kuznetsov and Mohri,
2016]. As the data generation process transitions to a stationary state with
an identical distribution at each step, this deviation $\Delta$ diminishes to
zero.

######  Assumption 4.

We assume the discrepancy measure that quantifies the deviation between the
intermediate steps (i.e., $i=1,...,N-1$) and the final step (i.e., $i=N$) data
distribution as

$\Delta=\sup_{f(\bm{\theta})}\Big{|}\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\left[\text{loss}^{0-1}_{i}(\bm{\theta}_{i-1})|\mathcal{D}_{1}^{i-1}\right]-\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\left[\text{loss}^{0-1}_{N}(\bm{\theta}_{i-1})|\mathcal{D}_{1}^{N-1}\right]\Big{|},$

where the supremum is on any model,
$\mathcal{D}_{1}^{i-1}=\\{(X_{1},y_{1}),\ldots,(X_{i-1},y_{i-1})\\}$ is the
sequence of data points before the $i$-th iteration, and
$\text{loss}^{0-1}_{i}(\bm{\theta})=\mathbf{1}\\{y_{i}f_{i}(\bm{\theta})<0\\}$
is 0-1 loss.

We introduce the feature-label alignment (FLA) score, which measures how well
the representations produced by different algorithms align with the ground-
truth labels.

######  Definition 1.

FLA is defined as
$\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}$, where
$\mathbf{J}=[\text{vec}(\nabla_{\theta}f_{i}(\bm{\theta}_{0}))]_{i=1}^{N}$ is
the gradient of different temporal graph algorithms computed on each
individual training example and $\bm{\theta}_{0}$ is the untrained weight
parameters initialized from a Gaussian distribution.

FLA has appeared in the convergence and generalization analysis of over-
parameterized neural networks [Arora et al., 2019]. In practice, FLA
quantifies the amount of perturbation we need on $\bm{\theta}$ along the
direction of
$\mathbf{J}=[\text{vec}(\nabla_{\theta}f_{i}(\bm{\theta}_{0}))]_{i=1}^{N}$ to
minimize the logistic loss, which could be used to capture the expressiveness
of different TGL algorithms, i.e., the smaller the perturbation, the better
the expressiveness. Detailed discussion can be found in Appendix I.1.
Computing FLA requires $\mathcal{O}(N^{2}|\bm{\theta}|)$ time complexity to
compute $\mathbf{J}\mathbf{J}^{\top}$ and $\mathcal{O}(N^{3})$ time complexity
to compute matrix inverse, where $N$ is the number of training data and
$|\bm{\theta}|$ is the number of weight parameters. In over-parameterized
regime, we assume $|\bm{\theta}|>N$, and we can compute
$\mathbf{J}\mathbf{J}^{\top}$ on a sampled subset of training data rather than
the full training data to make sure this assumption holds.

The following theorem establishes the generalization capability of different
TGL algorithms.

######  Theorem 1.

Given any $\delta\in(0,1/e]$, FLA-related constant
$R=\mathcal{O}(\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}})$,
and number of training iterations $N$ (one training example per iteration),
there exists $m^{\star}=\mathcal{O}(N^{2}/L^{2})\log(1/\delta)$ such that, if
hidden dimension $m\geq m^{\star}$ and using learning rate
$\eta=\frac{R}{mC\sqrt{2N}}$, with probability at least $1-\delta$ over the
randomness of $\bm{\theta}_{0}$ initialization, we can upper bound the
expected 0-1 error by

$\mathbb{E}[\text{loss}_{N}^{0-1}(\widetilde{\bm{\theta}})\leavevmode\nobreak\
|\leavevmode\nobreak\
\mathcal{D}_{1}^{N-1}]\leq\mathcal{O}\Big{(}\frac{DRC}{\sqrt{N}}\Big{)}+\mathcal{O}\Big{(}\sqrt{\frac{\ln(1/\delta)}{N}}\Big{)}+\Delta,$

where $\widetilde{\bm{\theta}}$ is uniformly sampled from
$\\{\bm{\theta}_{0},\ldots,\bm{\theta}_{N-1}\\}$,
$R=\mathcal{O}(\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}})$,
and the constant $C$ and $D$ of 1 the $L$-layer GNN are
$C=((1+3\rho)\tau)^{L-1}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
D=L$, 2 the $L$-step RNN are $C=(1+3\rho/\sqrt{2})^{L-1}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ D=L$, and 3 the memory-based method are
$C=\rho\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ D=4$.

In Theorem 1, we show that the generalization error decreases as the number of
training iterations $N$ increases, with a single data point being used for
training at each iteration. On the other hand, the generalization error
increases with respect to the number of layers/steps $L$ in GNN-/RNN-based
methods, the feature-label alignment constant $R$, the maximum Lipschitz
constant of the activation function $\rho$, and the graph convolution-related
constant $\tau$. In the following section, we delve into more details on the
effect of the number of layers/steps $L$ and the feature-label alignment
constant $R$. The proofs are deferred to Appendices C, D, and E, respectively.

###  3.3 Discussion on the generalization bound: insights and limitations

#### Dependency on depth and steps.

The generalization error of GNN- and RNN-based methods tends to increase as
the number of layers/steps $L$ increases. This partially explains why the
hyper-parameter selection on $L$ is usually small for those methods. For
example, _GNN-based method_ TGAT uses 2-layer GNN (i.e., $L=3$), _RNN-based
method_ CAW selects 3-steps RNN (i.e., $L=4$), and _GNN &RNN-based method_
DySAT uses 2-layer GNN and 3-steps RNN to achieve outstanding performance. On
the other hand, _memory-based method_ JODIE alleviates the dependency issue by
using memory blocks and can leverage all historical interactions for
prediction, which enjoys a generalization error independent of the number of
steps. However, since gradients cannot flow through the memory blocks due to
“stop gradient”, its expressive power may be lower than that of other methods,
which will be further elaborated when discussing the impact of feature-label
alignment. _Memory &GNN-based methods_ TGN and APAN alleviate the lack of
expressive power issue by applying a single-layer GNN on top of the memory
blocks.

#### Dependency on feature-label alignment (FLA).

Although the dependency on the number of layers/steps of GNN/RNN can partially
explain the performance disparity between these methods, it is still not clear
if using “stop gradient” in the memory-based method and the selection of input
data (e.g., using recent or uniformly sampled neighbors in a temporal graph)
can affect the model performance. In the following, we take a closer look at
the FLA score, which is inversely correlated to the generalization ability of
the TGL models, i.e., the smaller the FLA, the better the generalization
ability. According to Figure 3, we observe that 1 JODIE (memory-based) has a
relatively larger FLA score than most of the other TGL methods. This is
potentially due to “stop gradient” operation that prevents gradients from
flowing through the memory blocks and could potentially hurt the expressive
power. 2 APAN and TGN (memory&GNN-based) alleviate the expressive power
degradation issue by applying a single layer GNN on top of the memory blocks,
resulting in a smaller FLA than the pure memory-based method. 3 TGAT (GNN-
based) has a relatively smaller FLA score than other methods, which is
expected since GNN has achieved outstanding performance on static graphs. 4
DySAT (GNN&RNN-based) is originally designed for snapshot graphs, so its FLA
score might be highly dependent on the selection of time-span size when
converting a temporal graph to snapshot graphs. A non-optimal choice of time-
span might cause information loss, which partially explains why DySAT’s FLA is
large. Additionally, the selection of the input data also affects the FLA. We
will discuss further in the experimental section with empirical validation and
detailed analysis.

Figure 3: Comparison of FLA (y-axis) of different methods (x-axis) on real-
world datasets.

##  4 A simplified algorithm

Informed by our theoretical analysis, we introduce Simplified-Temporal-Graph-
Network (SToNe) that not only enjoys a small generalization error but also
empirically works well. The design of SToNe is guided by the following key
insights that are presented in Section 3.3:

  * •

Shallow and non-recursive network. This is because the generalization error
increases with respect to the number of layers/steps in GNN-/RNN-based
methods, which motivates us to consider a shallow and non-recursive neural
architecture to alleviate such dependency.

  * •

Selecting proper input data instead of using memory blocks. Although memory
blocks could alleviate the dependency of generalization error on the number of
layers/steps, it will also affect the FLA and hurt the generalization ability
of the models. As an alternative, we propose to capture the important
historical interactions by empirically selecting the proper input data.

To this end, we introduce the data preparation and the neural architecture in
Section 4.1, then highlight the key features of SToNe that can differentiate
itself from existing methods in Section 4.2.

###  4.1 Simplified temporal graph network: input data and neural architecture

Input data preparation. To compute the representation of node $v_{i}$ at time
$t$, we first identify the $K$ most recent nodes that have interacted with
$v_{i}$ prior to time $t$ and denote them as temporal neighbors
$\mathcal{N}_{K}^{t}(v_{i})$. Then, we sort all nodes inside
$\mathcal{N}_{K}^{t}(v_{i})$ by the descending temporal order. If a node
$v_{j}$ interacts with node $v_{i}$ multiple times, each interaction is
treated as a separate temporal neighbor. For example in Figure 2, for any time
$t>t_{6}$ and large enough constant $K$, we have
$\mathcal{N}_{K}^{t}(v_{4})=\\{(v_{5},t_{6}),(v_{2},t_{4}),(v_{2},t_{3})\\}$
in the descending temporal order. For each temporal neighbor
$v_{j}\in\mathcal{N}_{K}(v_{i},t)$, we represent its interaction with the
target node $v_{i}$ at time $t^{\prime}$ using a combination of edge features
$\mathbf{e}_{ij}(t^{\prime})\in\mathbb{R}^{d_{\text{e}}}$, time-encoding
$\bm{\psi}(t-t^{\prime})\in\mathbb{R}^{d_{\text{t}}}$, and node features
$\mathbf{x}_{i},\leavevmode\nobreak\
\mathbf{x}_{j}\in\mathbb{R}^{d_{\text{n}}}$, which is denoted as
$\mathbf{u}_{ij}(t^{\prime})=[\mathbf{e}_{ij}(t^{\prime})\leavevmode\nobreak\
||\leavevmode\nobreak\ \bm{\psi}(t-t^{\prime})\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{i}\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{j}]$. Then, we define
$\mathbb{H}_{i}(t)=\\{\mathbf{u}_{ij}(t^{\prime})\leavevmode\nobreak\
|\leavevmode\nobreak\ (v_{j},t^{\prime})\in\mathcal{N}_{K}^{t}(v_{i})\\}$ as
the set of all features that represent the temporal neighbors of node $v_{i}$
at time $t$, where all vectors in $\mathbb{H}_{i}(t)$ are ordered by the
descending temporal order. For example, the interactions between the temporal
neighbors of node $v_{4}$ at time $t>t_{6}$ to its root node in Figure 2 are
$\mathbb{H}_{4}(t)=\\{\mathbf{u}_{4,5}(t_{6}),\mathbf{u}_{4,2}(t_{4}),\mathbf{u}_{4,2}(t_{3})\\}$,
where $\mathbf{u}_{4,5}(t_{6})=[\mathbf{e}_{4,5}(t_{6})\leavevmode\nobreak\
||\leavevmode\nobreak\ \bm{\psi}(t-t_{6})\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{4}\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{5}],\leavevmode\nobreak\
\mathbf{u}_{4,2}(t_{4})=[\mathbf{e}_{4,2}(t_{4})\leavevmode\nobreak\
||\leavevmode\nobreak\ \bm{\psi}(t-t_{4})\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{4}\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{2}]$, and
$\mathbf{u}_{4,2}(t_{3})=[\mathbf{e}_{4,2}(t_{3})\leavevmode\nobreak\
||\leavevmode\nobreak\ \bm{\psi}(t-t_{3})\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{4}\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{x}_{2}]$. The time-encoding function in SToNe
is defined as $\bm{\psi}(t-t^{\prime})=\cos((t-t^{\prime})\mathbf{w})$, where
$\mathbf{w}=[\alpha^{-(i-1)/\beta}]_{i=1}^{d_{\text{t}}}$ is a fixed
$d_{\text{t}}$-dimensional vector and $\alpha=\beta=\sqrt{d_{\text{t}}}$.
Notice that a similar time-encoding function is used in other temporal graph
learning methods, e.g., [Kumar et al., 2019, Xu et al., 2020a, Rossi et al.,
2020, Wang et al., 2021e, Cong et al., 2023].

#### Encoding features via GNN.

SToNe is a GNN-based method with trainable parameters
$\bm{\theta}=\\{\bm{\alpha},\mathbf{W}^{(1)},\mathbf{W}^{(2)}\\}$, where
$\bm{\alpha}\in\mathbb{R}^{K}$,
$\mathbf{W}^{(1)}\in\mathbb{R}^{d_{\text{hid}}\times d_{\text{in}}}$, and
$\mathbf{W}^{(2)}\in\mathbb{R}^{d_{\text{out}}\times d_{\text{hid}}}$. Here
$K$ is the maximum temporal neighbor size we consider,
$d_{\text{in}}=d_{\text{e}}+d_{\text{t}}+2d_{\text{n}}$, and
$d_{\text{hid}},d_{\text{out}}$ are the dimensions of hidden and output
representations. The representation of node $v_{i}$ at time $t$ is computed by

$\mathbf{h}_{i}(t)=\mathbf{W}^{(2)}\text{LayerNorm}(\mathbf{z}_{i}(t)),\mathbf{z}_{i}(t)=\sigma\Big{(}\sum\nolimits_{k=1}^{|\mathbb{H}_{i}(t)|}\alpha_{k}\mathbf{W}^{(1)}[\mathbb{H}_{i}(t)]_{k}\Big{)}+\sum\nolimits_{k=1}^{|\mathbb{H}_{i}(t)|}[\mathbb{H}_{i}(t)]_{k},$
(1)

If the temporal neighbor size of node $v_{i}$ is less than $K$, i.e.,
$|\mathbb{H}_{i}(t)|<K$, we only need to update the first
$|\mathbb{H}_{i}(t)|$ entries of the vector $\bm{\alpha}$. Notice that the
number of parameters in SToNe is
$(K+d_{\text{hid}}d_{\text{in}}+d_{\text{hid}}d_{\text{out}})$, which is
usually fewer than other TGL algorithms given the same hidden dimension size.
We will report the number of parameters and its computational cost in the
experiment section.

Link prediction via MLP. When the downstream task is future link prediction,
we predict whether an interaction between nodes $v_{i},v_{j}$ occurs at time
$t$ by applying a 2-layer MLP model on $[\mathbf{h}_{i}(t)\leavevmode\nobreak\
||\leavevmode\nobreak\ \mathbf{h}_{j}(t)]\in\mathbb{R}^{2d_{\text{out}}}$. It
is worth noting that the same link classifier is used in almost all the
existing temporal graph learning methods [Kumar et al., 2019, Xu et al.,
2020a, Rossi et al., 2020, Wang et al., 2021e, Souza et al., 2022, Zhou et
al., 2022, Cong et al., 2023], including SToNe.

###  4.2 Comparison to existing methods

#### Comparison to TGAT.

GNN-based method TGAT uses 2-hop uniformly sampled neighbors and aggregates
the information using a 2-layer GAT [Veličković et al., 2017]. The neighbor
aggregation weights in TGAT are estimated by self-attention. In contrast,
SToNe uses 1-hop most recent neighbors and directly learns the neighbor
aggregation weights $\bm{\alpha}$ as shown in Eq. 1. Moreover, self-attention
in TGAT can be thought of as weighted average of neighbor information, while
SToNe can be thought of as sum aggregation, which can better distinguish
different temporal neighbors, and it is especially helpful when node and edge
features are lacking.

#### Comparison to TGN.

Memory&GNN-based method TGN uses 1-hop most recent temporal neighbors and
applies a self-attention module to the sampled temporal neighbors’ features
that are stored inside the memory blocks. In fact, SToNe can be thought of as
a special case of TGN, where we use the features in $\mathbb{H}_{i}(t)$
instead of the memory blocks and directly learn the neighbor aggregation
weight $\bm{\alpha}$ instead of using the self-attention aggregation as shown
in Eq. 1.

#### Comparison to GraphMixer.

SToNe could be think of as a simplified version of [Cong et al., 2023] that
addresses the high computation cost associated with the MLP-mixer used for
temporal aggregation [Tolstikhin et al., 2021]. Instead of relying on the MLP-
mixer, we introduce the aggregation vector $\bm{\alpha}$ and employ linear
functions parameterized by $\smash{\mathbf{W}^{(1)}}$ and
$\smash{\mathbf{W}^{(2)}}$ for aggregation. Additionally, we do not explicitly
model the graph structure through [Cong et al., 2023]’s node-encoder. Instead,
we implicitly capture the node features within the temporal interactions
present in $\mathbb{H}(t)$. Our experiments demonstrate that SToNe
significantly reduces the model complexity (Table 2) while achieving
comparable performance (Figure 4 and Table 1).

#### Temporal graph construction.

Most of the TGL methods [Kumar et al., 2019, Xu et al., 2020a, Rossi et al.,
2020, Sankar et al., 2020, Wang et al., 2021e] implements temporal graphs as
directed graph data structure with information only flowing from source to
destination nodes. However, we consider the temporal graph as an bi-directed
graph data structure by assuming that information also flow from destination
to source nodes for SToNe. This means that the “most recent 1-hop neighbors”
sampled for the two nodes on the “bi-directed” temporal graph should be
similar if two nodes are frequently connected in recent timestamps. Such
similarity provides information on whether two nodes are frequently connected
in recent timestamps, which is essential for temporal graph link prediction
[Cong et al., 2023]. For example, let us assume nodes $v_{i},v_{j}$ interacts
at time $t_{1},t_{2}$ in the temporal order. Then, given any timestamp
$t>t_{2}$ and a large enough $K$, the 1-hop temporal neighbors on the bi-
directed graph is
$\mathcal{N}_{K}^{t}(v_{i})=\\{(v_{j},t_{1}),(v_{j},t_{2})\\}$ and
$\mathcal{N}_{K}^{t}(v_{j})=\\{(v_{i},t_{1}),(v_{i},t_{2})\\}$, while on
directed graph is
$\mathcal{N}_{K}^{t}(v_{i})=\\{(v_{j},t_{1}),(v_{j},t_{2})\\}$ but
$\mathcal{N}_{K}^{t}(v_{j})=\emptyset$. Intuitively, if two nodes are
frequently connected in recent timestamps, they are also likely to be
connected in the near future. In the experiment section, we show that changing
the temporal graph from bi-directed to directed can negatively impact the
feature-label alignment and model performance.

##  5 Experiments

We compare SToNe with several TGL algorithms under the transductive learning
setting. We conduct experiments on 6 real-world datasets, i.e., Reddit, Wiki,
MOOC, LastFM, GDELT, and UCI. Similar to many existing works, we also use the
$70\%/15\%/15\%$ chronological splits for the train/validation/test sets. We
re-run the official released implementations on benchmark datasets and repeat
6 times with different random seeds. Please refer to Appendix A for experiment
setup details.

###  5.1 Experiment results

Comparison on average precision. We compare the average precision score with
baseline methods in Table 1. We observe that SToNe could achieve compatible or
even better performance on most datasets. In particular, the performance of
SToNe outperforms most of the baselines with a large margin on the LastFM and
UCI dataset. This is potentially because these two datasets lack node/edge
features and have larger average time-gap (see dataset statistics in Table 3).
Since baseline methods rely on RNN or self-attention to process the node/edge
features, they implicitly assume that the features exists and are “smooth” at
adjacent timestamps, which could generalize poorly when the assumptions are
violated. PINT addresses this issue by pre-processing the dataset and
generating its own positional encoding as the augment features, therefore it
achieves better performance than SToNe on UCI dataset. However, computing this
positional encoding is time-consuming and does not always perform well on
other datasets. For example, PINT requires more than $400$ hours to compute
the positional features on GDELT, which is infeasible. Additionally, the
performance of SToNe closely matches that of GraphMixer, but with
significantly lower model complexity, which will be demonstrated in Table 2.
SToNe exhibits improved and more consistent performance on UCI dataset,
because SToNe is less susceptible to overfitting on small dataset.

Table 1: Comparison on the average precision score for link prediction. We
highlight the best score in red, the second best score in blue, and the third
best in green.

| Reddit | Wiki | MOOC | LastFM | GDELT | UCI
---|---|---|---|---|---|---
JODIE | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.75}\pm 0.02$ | $98.94\pm 0.06$ | $98.99\pm 0.04$ | $79.41\pm 3.68$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}98.48}\pm 0.01$ | $73.93\pm 2.77$
TGAT | $99.56\pm 0.04$ | $98.69\pm 0.10$ | $99.28\pm 0.07$ | $75.16\pm 0.10$ | $96.46\pm 0.04$ | $82.49\pm 0.67$
TGSRec | $95.21\pm 0.08$ | $91.64\pm 0.12$ | $83.62\pm 0.34$ | $76.91\pm 0.87$ | $97.03\pm 0.61$ | $76.64\pm 0.54$
TGN | $98.83\pm 0.01$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.61}\pm 0.05$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.63}\pm 0.06$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}91.04}\pm 2.18$ | $98.33\pm 0.09$ | $79.91\pm 1.59$
APAN | $99.36\pm 0.17$ | $98.99\pm 0.14$ | $99.02\pm 0.11$ | $73.18\pm 7.72$ | $98.01\pm 0.26$ | $61.73\pm 5.96$
PINT | $99.03\pm 0.01$ | $98.78\pm 0.10$ | $87.24\pm 0.73$ | $88.06\pm 0.71$ | - | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}96.22}\pm 0.08$
CAW-attn | $98.51\pm 0.02$ | $97.95\pm 0.03$ | $63.07\pm 0.82$ | $76.31\pm 0.10$ | $95.06\pm 0.11$ | $92.16\pm 0.19$
DySAT | $98.55\pm 0.01$ | $96.64\pm 0.03$ | $98.76\pm 0.11$ | $76.28\pm 0.04$ | $98.17\pm 0.01$ | $80.43\pm 0.36$
GraphMixer | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.93}\pm 0.01$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.85}\pm 0.01$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.91}\pm 0.01$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}96.31}\pm 0.02$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}98.89}\pm 0.02$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}92.39}\pm 2.15$
SToNe | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.89}\pm 0.00$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.85}\pm 0.05$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.88}\pm 0.04$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}95.74}\pm 0.13$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.11}\pm 0.03$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}94.60}\pm 0.31$

Comparison on generalization performance. To validate the generalization
performance of SToNe, we compare the average precision score and the
generalization gap in Figure 4. The generalization gap is defined as the
absolute difference between the training and validation average precision
scores. Our results show that SToNe, similar to GraphMixer, consistently
achieves a higher average precision score faster than other baselines, has a
smaller generalization gap, and has relatively stable performance across
different epochs. This suggests that SToNe has better generalization
performance compared to the baselines. In particular on the UCI dataset, SToNe
is less prone to overfit and its generalization gap increases more slowly than
all other baseline methods.

Figure 4: Comparison of the average prevision of validation set and
generalization gap of different methods on real-world datasets.

Comparison on model complexity. We compare the number of parameters (including
trainable weight parameters and memory block size) and wall-clock time per
epoch during the training phase in Table 2. Our results show that SToNe has
fewer parameters than all the baselines, and its computation time is also
faster than most baseline methods. Note that some baselines also require
significant time for data preprocessing, which is not included in Table 2. For
example, PINT takes more than $84$ hours to pre-compute the positional
encoding on the Reddit dataset, which is significantly longer than its per
epoch training time. In contrast, SToNe does not require computing augmented
features based on the temporal graph and therefore does not have this pre-
processing time.

Table 2:  Comparison of the number of parameters ($\times 10^{5}$ parameters)
and wall-clock time (second) of single epoch of training in the format “Number
of parameters (Wall-clock time)”.

| JODIE | TGAT | TGSRec | TGN | APAN | PINT | CAW | DySAT | GraphMixer | SToNe
---|---|---|---|---|---|---|---|---|---|---
Reddit | 11.6 (5s) | 2.1 (15s) | 5.1 (538s) | 14.1 (8s) | 12.3 (13s) | 17.1 (436s) | 43.3 (1,930s) | 4.3 (33s) | 23.3 (12s) | 0.58 (8s)
Wiki | 9.8 (2s) | 2.1 (4s) | 4.8 (157s) | 12.3 (2s) | 10.6 (4s) | 17.9 (93s) | 43.3 (282s) | 4.3 (6s) | 19.8 (3s) | 0.58 (2s)
MOOC | 7.8 (4s) | 1.4 (8s) | 3.9 (656s) | 10.0 (5s) | 8.4 (9s) | 10.1 (157s) | 43.3 (653s) | 2.3 (16d) | 15.3 (7s) | 0.41 (5s)
LastFM | 2.6 (11s) | 1.4 (28s) | 2.3 (1,810s) | 4.8 (15s) | 3.2 (28s) | 4.9 (440s) | 43.3 (1,832s) | 2.3 (41s) | 5.2 (21s) | 0.41 (16s)

More experiment results. Due to the space limit, we defer more experiment
results to Appendix B, such as performance evaluation under different metrics
(e.g., AUC, RecallK, and MRR), ablation study on the effect of neighbor
selection (e.g., recent vs randomly sampled neighbors), etc.

##  6 Conclusion

In this paper, we study the generalization ability of various TGL algorithms.
We reveal the relationship between the generalization error to “the number of
layers/steps” and “the feature-label alignment”. Guided by our analysis, we
propose SToNe. Extensive experiments are conducted on real-world datasets to
show that SToNe enjoys a smaller generalization error, better performance, and
lower complexity. These results provide a deeper insight into TGL from a
theoretical perspective and are beneficial to design practically efficacious
and simpler TGL algorithms for future studies.

## Acknowledgements

This work is supported in part by NSF Award #2008398, #2134079, and #1939725.

## References

  * [Abboud et al., 2020] Abboud, R., Ceylan, I. I., Grohe, M., and Lukasiewicz, T. (2020).  The surprising power of graph neural networks with random node initialization.  arXiv preprint arXiv:2010.01179.
  * [Arora et al., 2019] Arora, S., Du, S., Hu, W., Li, Z., and Wang, R. (2019).  Fine-grained analysis of optimization and generalization for overparameterized two-layer neural networks.  In International Conference on Machine Learning.
  * [Arora et al., 2022] Arora, S., Li, Z., and Panigrahi, A. (2022).  Understanding gradient descent on the edge of stability in deep learning.  In International Conference on Machine Learning, pages 948–1024. PMLR.
  * [Bouritsas et al., 2022] Bouritsas, G., Frasca, F., Zafeiriou, S. P., and Bronstein, M. (2022).  Improving graph neural network expressivity via subgraph isomorphism counting.  IEEE Transactions on Pattern Analysis and Machine Intelligence.
  * [Cai et al., 2022] Cai, B., Xiang, Y., Gao, L., Zhang, H., Li, Y., and Li, J. (2022).  Temporal knowledge graph completion: A survey.  arXiv preprint arXiv:2201.08236.
  * [Cao and Gu, 2019] Cao, Y. and Gu, Q. (2019).  Generalization bounds of stochastic gradient descent for wide and deep neural networks.  Advances in neural information processing systems.
  * [Cesa-Bianchi et al., 2004] Cesa-Bianchi, N., Conconi, A., and Gentile, C. (2004).  On the generalization ability of on-line learning algorithms.  IEEE Transactions on Information Theory, 50(9):2050–2057.
  * [Chi et al., 2022] Chi, H., Xu, H., Fu, H., Liu, M., Zhang, M., Yang, Y., Hao, Q., and Wu, W. (2022).  Long short-term preference modeling for continuous-time sequential recommendation.  arXiv preprint arXiv:2208.00593.
  * [Cong et al., 2021] Cong, W., Ramezani, M., and Mahdavi, M. (2021).  On provable benefits of depth in training graph convolutional networks.  In Advances in Neural Information Processing Systems.
  * [Cong et al., 2023] Cong, W., Zhang, S., Kang, J., Yuan, B., Wu, H., Zhou, X., Tong, H., and Mahdavi, M. (2023).  Do we really need complicated model architectures for temporal networks?  In The Eleventh International Conference on Learning Representations.
  * [Du et al., 2019] Du, S. S., Hou, K., Salakhutdinov, R., Póczos, B., Wang, R., and Xu, K. (2019).  Graph neural tangent kernel: Fusing graph neural networks with graph kernels.  In Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada.
  * [Du et al., 2018] Du, S. S., Zhai, X., Poczos, B., and Singh, A. (2018).  Gradient descent provably optimizes over-parameterized neural networks.  arXiv preprint arXiv:1810.02054.
  * [Fan et al., 2021] Fan, Z., Liu, Z., Zhang, J., Xiong, Y., Zheng, L., and Yu, P. S. (2021).  Continuous-time sequential recommendation with temporal graph collaborative transformer.  In Proceedings of the 30th ACM International Conference on Information & Knowledge Management.
  * [Fey and Lenssen, 2019] Fey, M. and Lenssen, J. E. (2019).  Fast graph representation learning with PyTorch Geometric.  In ICLR Workshop on Representation Learning on Graphs and Manifolds.
  * [Gao and Ribeiro, 2022] Gao, J. and Ribeiro, B. (2022).  On the equivalence between temporal and static equivariant graph representations.  In International Conference on Machine Learning.
  * [Gao et al., 2023] Gao, J., Zhou, Y., and Ribeiro, B. (2023).  Double permutation equivariance for knowledge graph completion.  arXiv preprint arXiv:2302.01313.
  * [Garg et al., 2020] Garg, V., Jegelka, S., and Jaakkola, T. (2020).  Generalization and representational limits of graph neural networks.  In International Conference on Machine Learning.
  * [Goyal et al., 2018] Goyal, P., Chhetri, S. R., and Canedo, A. (2018).  dyngraph2vec: Capturing network dynamics using dynamic graph representation learning.  arXiv preprint arXiv:1809.02657.
  * [Hajiramezanali et al., 2019] Hajiramezanali, E., Hasanzadeh, A., Narayanan, K., Duffield, N., Zhou, M., and Qian, X. (2019).  Variational graph recurrent neural networks.  Advances in neural information processing systems.
  * [Hron et al., 2020] Hron, J., Bahri, Y., Sohl-Dickstein, J., and Novak, R. (2020).  Infinite attention: Nngp and ntk for deep attention networks.  In International Conference on Machine Learning, pages 4376–4386. PMLR.
  * [Hu et al., 2020] Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., and Leskovec, J. (2020).  Open graph benchmark: Datasets for machine learning on graphs.  Advances in neural information processing systems.
  * [Huang et al., 2023] Huang, S., Poursafaei, F., Danovitch, J., Fey, M., Hu, W., Rossi, E., Leskovec, J., Bronstein, M. M., Rabusseau, G., and Rabbany, R. (2023).  Temporal graph benchmark for machine learning on temporal graphs.  In Thirty-seventh Conference on Neural Information Processing Systems Datasets and Benchmarks Track.
  * [Kumar et al., 2019] Kumar, S., Zhang, X., and Leskovec, J. (2019).  Predicting dynamic embedding trajectory in temporal interaction networks.  In Proceedings of the 25th ACM SIGKDD international conference on Knowledge discovery and data mining.
  * [Kuznetsov and Mohri, 2015] Kuznetsov, V. and Mohri, M. (2015).  Learning theory and algorithms for forecasting non-stationary time series.  Advances in neural information processing systems, 28.
  * [Kuznetsov and Mohri, 2016] Kuznetsov, V. and Mohri, M. (2016).  Time series prediction and online learning.  In Conference on Learning Theory, pages 1190–1213. PMLR.
  * [Leblay and Chekol, 2018] Leblay, J. and Chekol, M. W. (2018).  Deriving validity time in knowledge graph.  In Companion Proceedings of the The Web Conference 2018.
  * [Li et al., 2020] Li, P., Wang, Y., Wang, H., and Leskovec, J. (2020).  Distance encoding: Design provably more powerful neural networks for graph representation learning.  Advances in Neural Information Processing Systems.
  * [Liao et al., 2020] Liao, R., Urtasun, R., and Zemel, R. (2020).  A pac-bayesian approach to generalization bounds for graph neural networks.  arXiv preprint arXiv:2012.07690.
  * [Maskey et al., 2022] Maskey, S., Lee, Y., Levie, R., and Kutyniok, G. (2022).  Stability and generalization capabilities of message passing graph neural networks.  arXiv preprint arXiv:2202.00645.
  * [Nguyen, 2021] Nguyen, Q. (2021).  On the proof of global convergence of gradient descent for deep relu networks with linear widths.  In International Conference on Machine Learning.
  * [Oono and Suzuki, 2020] Oono, K. and Suzuki, T. (2020).  Optimization and generalization analysis of transduction through gradient boosting and application to multi-scale graph neural networks.  Advances in Neural Information Processing Systems.
  * [Pareja et al., 2020] Pareja, A., Domeniconi, G., Chen, J., Ma, T., Suzumura, T., Kanezashi, H., Kaler, T., Schardl, T., and Leiserson, C. (2020).  Evolvegcn: Evolving graph convolutional networks for dynamic graphs.  In Proceedings of the AAAI Conference on Artificial Intelligence.
  * [Rossi et al., 2020] Rossi, E., Chamberlain, B., Frasca, F., Eynard, D., Monti, F., and Bronstein, M. (2020).  Temporal graph networks for deep learning on dynamic graphs.  In ICML 2020 Workshop on Graph Representation Learning.
  * [Sankar et al., 2020] Sankar, A., Wu, Y., Gou, L., Zhang, W., and Yang, H. (2020).  Dysat: Deep neural representation learning on dynamic graphs via self-attention networks.  In Proceedings of the 13th International Conference on Web Search and Data Mining.
  * [Shi et al., 2020] Shi, Y., Huang, Z., Feng, S., Zhong, H., Wang, W., and Sun, Y. (2020).  Masked label prediction: Unified message passing model for semi-supervised classification.  arXiv preprint arXiv:2009.03509.
  * [Souza et al., 2022] Souza, A. H., Mesquita, D., Kaski, S., and Garg, V. (2022).  Provably expressive temporal graph networks.  arXiv preprint arXiv:2209.15059.
  * [Srinivasan and Ribeiro, 2019] Srinivasan, B. and Ribeiro, B. (2019).  On the equivalence between positional node embeddings and structural graph representations.  arXiv preprint arXiv:1910.00452.
  * [Tian et al., 2021] Tian, S., Wu, R., Shi, L., Zhu, L., and Xiong, T. (2021).  Self-supervised representation learning on dynamic graphs.  CIKM ’21.
  * [Tolstikhin et al., 2021] Tolstikhin, I. O., Houlsby, N., Kolesnikov, A., Beyer, L., Zhai, X., Unterthiner, T., Yung, J., Steiner, A., Keysers, D., Uszkoreit, J., et al. (2021).  Mlp-mixer: An all-mlp architecture for vision.  Advances in neural information processing systems, 34:24261–24272.
  * [Trivedi et al., 2019] Trivedi, R., Farajtabar, M., Biswal, P., and Zha, H. (2019).  Dyrep: Learning representations over dynamic graphs.  In International conference on learning representations.
  * [Veličković et al., 2017] Veličković, P., Cucurull, G., Casanova, A., Romero, A., Lio, P., and Bengio, Y. (2017).  Graph attention networks.  arXiv preprint arXiv:1710.10903.
  * [Verma and Zhang, 2019] Verma, S. and Zhang, Z. (2019).  Stability and generalization of graph convolutional neural networks.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD 2019, Anchorage, AK, USA, August 4-8, 2019.
  * [Vershynin, 2018] Vershynin, R. (2018).  High-dimensional probability: An introduction with applications in data science.
  * [Wang et al., 2022] Wang, H., Yin, H., Zhang, M., and Li, P. (2022).  Equivariant and stable positional encoding for more powerful graph neural networks.  arXiv preprint arXiv:2203.00199.
  * [Wang et al., 2021a] Wang, L., Chang, X., Li, S., Chu, Y., Li, H., Zhang, W., He, X., Song, L., Zhou, J., and Yang, H. (2021a).  Tcl: Transformer-based dynamic graph modelling via contrastive learning.  arXiv preprint arXiv:2105.07944.
  * [Wang et al., 2021b] Wang, X., Lyu, D., Li, M., Xia, Y., Yang, Q., Wang, X., Wang, X., Cui, P., Yang, Y., Sun, B., et al. (2021b).  Apan: Asynchronous propagation attention network for real-time temporal graph embedding.  In Proceedings of the 2021 International Conference on Management of Data.
  * [Wang et al., 2021c] Wang, Y., Cai, Y., Liang, Y., Ding, H., Wang, C., Bhatia, S., and Hooi, B. (2021c).  Adaptive data augmentation on temporal graphs.  Advances in Neural Information Processing Systems.
  * [Wang et al., 2021d] Wang, Y., Cai, Y., Liang, Y., Ding, H., Wang, C., and Hooi, B. (2021d).  Time-aware neighbor sampling for temporal graph networks.  arXiv preprint arXiv:2112.09845.
  * [Wang et al., 2021e] Wang, Y., Chang, Y.-Y., Liu, Y., Leskovec, J., and Li, P. (2021e).  Inductive representation learning in temporal networks via causal anonymous walks.  arXiv preprint arXiv:2101.05974.
  * [Xu et al., 2020a] Xu, D., Ruan, C., Korpeoglu, E., Kumar, S., and Achan, K. (2020a).  Inductive representation learning on temporal graphs.  arXiv preprint arXiv:2002.07962.
  * [Xu et al., 2018] Xu, K., Li, C., Tian, Y., Sonobe, T., Kawarabayashi, K.-i., and Jegelka, S. (2018).  Representation learning on graphs with jumping knowledge networks.  In International conference on machine learning.
  * [Xu et al., 2020b] Xu, K., Li, J., Zhang, M., Du, S. S., ichi Kawarabayashi, K., and Jegelka, S. (2020b).  What can neural networks reason about?  In International Conference on Learning Representations.
  * [Xu et al., 2020c] Xu, K., Li, J., Zhang, M., Du, S. S., Kawarabayashi, K., and Jegelka, S. (2020c).  What can neural networks reason about?  In 8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020.
  * [Xu et al., 2021a] Xu, K., Zhang, M., Jegelka, S., and Kawaguchi, K. (2021a).  Optimization of graph neural networks: Implicit acceleration by skip connections and more depth.  In International Conference on Machine Learning, pages 11592–11602. PMLR.
  * [Xu et al., 2021b] Xu, K., Zhang, M., Li, J., Du, S. S., Kawarabayashi, K.-I., and Jegelka, S. (2021b).  How neural networks extrapolate: From feedforward to graph neural networks.  In International Conference on Learning Representations.
  * [Yuan and Li, 2021] Yuan, H. and Li, G. (2021).  A survey of traffic prediction: from spatio-temporal data to intelligent transportation.  Data Science and Engineering.
  * [Zhang et al., 2021] Zhang, Q., Yu, K., Guo, Z., Garg, S., Rodrigues, J. J., Hassan, M. M., and Guizani, M. (2021).  Graph neural network-driven traffic forecasting for the connected internet of vehicles.  IEEE Transactions on Network Science and Engineering.
  * [Zhou et al., 2022] Zhou, H., Zheng, D., Nisa, I., Ioannidis, V., Song, X., and Karypis, G. (2022).  Tgl: A general framework for temporal gnn training on billion-scale graphs.  arXiv preprint arXiv:2203.14883.
  * [Zhou and Wang, 2021] Zhou, X. and Wang, H. (2021).  The generalization error of graph convolutional networks may enlarge with more layers.  Neurocomputing.
  * [Zhu et al., 2022] Zhu, Z., Liu, F., Chrysos, G. G., and Cevher, V. (2022).  Generalization properties of nas under activation and skip connection search.  arXiv preprint arXiv:2209.07238.

Organization. The supplementary material is organized as follows:

  * •

In Section A, we provide an overview on the hardware and software that used in
the experiments, including details on the datasets, as well as the
implementation of the baseline algorithms and SToNe.

  * •

In Section B, we present additional experimental results using different
evaluation metrics (AUC, Recall@K, and MRR), as well as ablation studies to
further analyze the contributions of different components in SToNe.

  * •

In Section C, Section D, and Section E, we provide the thorough proof on the
generalization analysis of GNN-based, RNN-based, and memory-based TGL
algorithms.

  * •

In Section H, we provide a summary of additional related works on the
expressive power and generalization of graph learning algorithms, as well as
additional TGL algorithms.

  * •

In Section I, we provide further discussions and clarifications on several
important details of this paper, e.g., details on related works and
backgrounds, discussions on our theoretical limitations, and future
clarifications on experiment settings and results.

##  Appendix A Experiment setup details

###  A.1 Hardware specification and environment

Our experiments are conducted on a single machine with an Intel i9-10850K
processor, an Nvidia RTX 3090 GPU, and 64GB of RAM. The experiments are
implemented in Python 3.8 using PyTorch 1.12.1 on CUDA 11.6 and the temporal
graph learning framework [Zhou et al., 2022]. To ensure a fair comparison, all
experiments are repeated 6 times with random seeds $\\{0,1,2,3,4,5\\}$.

###  A.2 Details on datasets

We conduct experiments on the following 6 datasets, where the detailed dataset
statistics are summarized in Table 3. The download links for all datasets can
be found in the code repository.

  * •

Reddit dataset consists of one month of posts made by users on subreddits. The
link feature is extracted by converting the text of each post into a feature
vector. Reddit dataset has been previously used in existing works such as
[Kumar et al., 2019, Xu et al., 2020a, Rossi et al., 2020, Souza et al., 2022,
Wang et al., 2021e, Zhou et al., 2022].

  * •

Wiki dataset consists of one month of edits made by edits on Wikipedia pages.
The link feature is extracted by converting the edit test into an LIWC-feature
vector. Wiki dataset has been previously used in existing works such as [Kumar
et al., 2019, Xu et al., 2020a, Rossi et al., 2020, Souza et al., 2022, Wang
et al., 2021e, Zhou et al., 2022].

  * •

LastFM dataset consists of one month of who-listen-to-which song information.
LastFM dataset has been previously used in existing works such as [Kumar et
al., 2019, Souza et al., 2022, Zhou et al., 2022].

  * •

MOOC dataset consists of actions done by students on a MOOC online course.
MOOC dataset has been previously used in existing works such as [Kumar et al.,
2019, Souza et al., 2022, Zhou et al., 2022].

  * •

GDELT dataset is a temporal knowledge graph dataset that originates from the
Event Database. The Event Database records events from around the world as
reported in news articles. It has been previously used in works such as [Zhou
et al., 2022, Cong et al., 2023]. We follow [Cong et al., 2023] to subsample 1
temporal link per 100 continuous temporal link because the original dataset is
too big to fit into CPU RAM memory for single-machine training.

  * •

UCI dataset is a publicly available communication network dataset that
includes email interactions between core employees and messages sent between
peer users on an online social network platform. It has been previously used
in works such as [Sankar et al., 2020, Souza et al., 2022].

Table 3: Dataset statistic.

| $|\mathcal{V}|$ | $|\mathcal{E}|$ | Avg time-gap | $\text{dim}(\mathbf{x}_{i}^{\text{node}})$ | $\text{dim}(\mathbf{x}_{ij}^{\text{link}})$ | Node | Link | Time
---|---|---|---|---|---|---|---|---
Reddit | 10,984 | 672,447 | 4 | 0 | 172 |  |  |
Wiki | 9,227 | 157,474 | 17 | 0 | 172 |  |  |
MOOC | 7,144 | 411,749 | 3.6 | 0 | 0 |  |  |
LastFM | 1,980 | 1,293,103 | 106 | 0 | 0 |  |  |
GDELT | 8,831 | 1,912,909 | 0.1 | 413 | 186 |  |  |
UCI | 1,900 | 59,835 | 279 | 0 | 0 |  |  |

###  A.3 Details on baseline implementations

The implementations of JODIE, DySAT, TGAT, TGN, and APAN are obtained from the
TGL framework [Zhou et al., 2022] at TGL-code111https://github.com/amazon-
research/tgl. This framework’s implementation has been found to achieve better
overall scores than the original implementations of these baselines.

The implementation of CAWs-attn is obtained from their official implementation
at CAW-code222https://github.com/snap-stanford/CAW.

The implementation of TGSRec is obtained from TGSRec-
code333https://github.com/DyGRec/TGSRec.

The implementation of PINT is obtained from PINT-
code444https://github.com/AaltoPML/PINT.

The implementation of GraphMixer is obtained from GraphMixer-
code555https://github.com/CongWeilin/GraphMixer.

We follow their instructions that are provided in the code repository for
hyper-parameter selection. We directly test with their official
implementations by changing our data structure to their required structure and
use their default hyper-parameters.

###  A.4 Details on SToNe implementations

We implement the proposed method SToNe under the TGL framework [Zhou et al.,
2022] and use the same hyper-parameters on all datasets (e.g., learning rate
$0.0001$, weight decay $10^{-6}$, batch size $600$, hidden dimension $100$) to
ensure a fair comparison. The models are trained until the validation error
did not decrease for 20 epochs. One of the most important dataset-dependent
hyper-parameter is the number of temporal graph neighbors $K$. The selection
of $K$ can be found in our repository.

##  Appendix B More experiment results

###  B.1 Transductive learning with AUC as evaluation metric

AUC (Under the ROC Curve) is one of the most widely accepted evaluation
metrics for link prediction, which has been used in many existing works [Xu et
al., 2020a, Rossi et al., 2020]. In the following, we compare the AUC score of
SToNe with baselines in Table 4. We observe that SToNe performs better than
the baselines on most datasets. While our performance is slightly lower than
PINT on the UCI dataset, PINT requires a significant amount of time to pre-
compute the positional features as augmented features (about 3 hours for the
UCI dataset) and the use of positional encoding in PINT does not always
perform well on other datasets. Besides, we do not report the results of PINT
on the GDELT dataset because it requires more than $400$ hours to compute the
positional features for training. Additionally, the performance of SToNe
closely matches that of GraphMixer, but with significantly lower model
complexity, which will be demonstrated in Table 2. SToNe exhibits improved and
more consistent performance compared on UCI dataset, which is attributed to
its small dataset size and SToNe is less susceptible to overfitting on small
dataset. Please also refer to our discussion on computational time and average
precision score in Section 5.1 for more details.

Table 4: Comparison on the AUC score for link prediction. We highlight the
best score in red, the second best score in blue, and the third best in green.

| Reddit | Wiki | MOOC | LastFM | GDELT | UCI
---|---|---|---|---|---|---
JODIE | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.78}\pm 0.01$ | $99.10\pm 0.05$ | $99.46\pm 0.07$ | $80.48\pm 5.84$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}98.66}\pm 0.01$ | $76.15\pm 3.28$
TGAT | $99.59\pm 0.04$ | $98.78\pm 0.08$ | $99.53\pm 0.03$ | $78.49\pm 0.09$ | $97.21\pm 0.02$ | $83.30\pm 0.40$
TGSRec | $94.74\pm 0.20$ | $91.32\pm 0.19$ | $80.70\pm 2.31$ | $76.66\pm 1.54$ | $96.72\pm 0.42$ | $81.00\pm 0.58$
TGN | $98.86\pm 0.01$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.62}\pm 0.05$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}99.77}\pm 0.04$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}92.54}\pm 2.24$ | $98.55\pm 0.05$ | $75.32\pm 2.34$
APAN | $99.50\pm 0.11$ | $99.12\pm 0.11$ | $99.43\pm 0.08$ | $78.15\pm 6.99$ | $98.39\pm 0.15$ | $64.39\pm 7.42$
PINT | $99.03\pm 0.01$ | $98.78\pm 0.10$ | $87.24\pm 0.73$ | $88.06\pm 0.71$ | - | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}94.22}\pm 0.58$
CAW-attn | $98.30\pm 0.01$ | $97.89\pm 0.02$ | $63.95\pm 0.81$ | $72.93\pm 0.54$ | $95.13\pm 0.11$ | $92.10\pm 0.10$
DySAT | $98.43\pm 0.01$ | $96.87\pm 0.03$ | $99.25\pm 0.05$ | $73.93\pm 0.08$ | $98.38\pm 0.01$ | $78.70\pm 0.05$
GraphMixer | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.94}\pm 0.01$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.82}\pm 0.01$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.93}\pm 0.01$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}97.38}\pm 0.02$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}98.89}\pm 0.02$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}91.74}\pm 2.69$
SToNe | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.91}\pm 0.00$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.85}\pm 0.05$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}99.91}\pm 0.02$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}97.16}\pm 0.06$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}99.26}\pm 0.02$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}93.14}\pm 0.67$

###  B.2 Ablation study on the effect of neighbor selection

We study the effect of model input selection on feature-label alignment and
average precision scores. As shown in Figure 5, changing the default setting
of “recent 1-hop neighbors” in SToNe to either “recent 2-hop neighbors” or
“random 1-hop neighbor” increases feature-label alignment score, which has a
inverse correlation to the generalization ability, and decreases average
precision scores.

Figure 5: Comparison of feature-label alignment (y-axis) and average precision
score (in red text) of different model input selection (x-axis) on real-world
datasets.

###  B.3 Directed versus bi-directed graph.

We study the impact of using directed/bi-directed graph data structure on the
feature-label alignment (FLA) and average precision score. As shown in Figure
6, changing from a bi-directed to a directed graph will increase the FLA score
and lead to a significant decrease in model performance. This is because “the
recent 1-hop neighbors” sampled for the two nodes on the bi-directed graph
could provide information on whether two nodes are frequently connected in the
last few timestamps, which is essential for temporal graph link prediction.
However, such type of information is missing if using directed graph as
discussed in the last paragraph of Section 4.2.

Figure 6: Compare the feature-label alignment (y-axis) and average precision
(red text) with directed/undirected temporal graph.

###  B.4 Transductive learning with Recall@K and MRR as evaluation metrics

Recall@K and MRR (mean reciprocal rank) are popular evaluation metrics
commonly used in the real-world recommendation system. Higher values indicate
better model performance. Our implementation of Recall@K and MRR follows the
implementation in [Hu et al., 2020, Cong et al., 2023]: we first sample 100
negative destination nodes for each source node of a temporal link node pair,
then our goal is to rank the positive temporal link node pairs higher than 100
negative destination nodes. As shown in Table 5, SToNe and GraphMixer has
higher Recall@K and MRR scores than other baselines, indicating that SToNe has
more confidence in the estimated categories. This might because SToNe has
better generalization ability and higher confidence in its predictions when
using the knowledge learned from the training set.

Table 5: Comparison on the Recall@1, Recall@5 and, MRR. For exch row, the best
score is in red text, second best score is in blue text, and the first best
score is in green.

|  | JODIE | DySAT | TGAT | TGN | GraphMixer | SToNe
---|---|---|---|---|---|---|---
Reddit | Recall@1 | $0.8773\pm 0.0121$ | $0.8454\pm 0.0007$ | $0.8428\pm 0.0130$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9294}\pm 0.0028$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9999}\pm 0.0000$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9998}\pm 0.0001$
Recall@5 | $0.9712\pm 0.0012$ | $0.9362\pm 0.0006$ | $0.9531\pm 0.0077$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9823}\pm 0.0012$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1.0000}\pm 0.0000$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1.0000}\pm 0.0000$
MRR | $0.8425\pm 0.0152$ | $0.8153\pm 0.0016$ | $0.8219\pm 0.0116$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9072}\pm 0.0040$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9965}\pm 0.0000$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9958}\pm 0.0002$
Wiki | Recall@1 | $0.8347\pm 0.0051$ | $0.7689\pm 0.0043$ | $0.7139\pm 0.0150$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.8526}\pm 0.0050$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9954}\pm 0.0002$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9955}\pm 0.0002$
Recall@5 | $0.9286\pm 0.0033$ | $0.8763\pm 0.0080$ | $0.8915\pm 0.0090$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9293}\pm 0.0026$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9971}\pm 0.0002$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9980}\pm 0.0001$
MRR | $0.8145\pm 0.0048$ | $0.7450\pm 0.0084$ | $0.6999\pm 0.0128$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.8343}\pm 0.0063$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9875}\pm 0.0004$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9827}\pm 0.0009$
MOOC | Recall@1 | $0.8637\pm 0.0130$ | $0.8337\pm 0.0059$ | $0.9066\pm 0.0072$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9314}\pm 0.0175$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9994\pm 0.0001}$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9994}\pm 0.0001$
Recall@5 | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9981}\pm 0.0002$ | $0.9981\pm 0.0002$ | $0.9927\pm 0.0013$ | $0.9932\pm 0.0009$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9999}\pm 0.0001$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9999}\pm 0.0000$
MRR | $0.7675\pm 0.0151$ | $0.7451\pm 0.0058$ | $0.8234\pm 0.0083$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.8602}\pm 0.0263$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9911}\pm 0.0004$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9910}\pm 0.0006$
LastFM | Recall@1 | $0.1684\pm 0.0589$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.2773}\pm 0.0035$ | $0.1770\pm 0.0041$ | $0.2746\pm 0.1279$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9940}\pm 0.0013$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9951}\pm 0.0002$
Recall@5 | $0.2982\pm 0.0766$ | $0.3886\pm 0.0050$ | $0.2617\pm 0.0043$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.4312}\pm 0.2001$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9997}\pm 0.0002$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}1.0000}\pm 0.0000$
MRR | $0.2071\pm 0.0536$ | $0.2892\pm 0.0074$ | $0.1998\pm 0.0042$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.2995}\pm 0.1278$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9620}\pm 0.0015$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9592}\pm 0.0011$
GDELT | Recall@1 | $0.7533\pm 0.0053$ | $0.3503\pm 0.0138$ | $0.0182\pm 0.0053$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.7595}\pm 0.0047$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9115}\pm 0.0028$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9117}\pm 0.0027$
Recall@5 | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.9333}\pm 0.0011$ | $0.8239\pm 0.0068$ | $0.0576\pm 0.0134$ | $0.9318\pm 0.0023$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.9891}\pm 0.0011$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9899}\pm 0.0009$
MRR | $0.7108\pm 0.0076$ | $0.4146\pm 0.0103$ | $0.0868\pm 0.0061$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.7302}\pm 0.0044$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.8622}\pm 0.0037$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.8627}\pm 0.0032$
UCI | Recall@1 | $0.5736\pm 0.0068$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.6039}\pm 0.0878$ | $0.1872\pm 0.0292$ | $0.4732\pm 0.0434$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.7989}\pm 0.1119$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.8995}\pm 0.0069$
Recall@5 | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.7814}\pm 0.0022$ | $0.6682\pm 0.0889$ | $0.3036\pm 0.0339$ | $0.6384\pm 0.0396$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.8571}\pm 0.0771$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.9435}\pm 0.0034$
MRR | $0.5572\pm 0.0087$ | ${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}0.5656}\pm 0.0816$ | $0.2296\pm 0.0253$ | $0.4799\pm 0.0354$ | ${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}0.7677}\pm 0.1093$ | ${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}0.8628}\pm 0.0087$

###  B.5 Ablation study on the effect of using self-attention in SToNe

In this section, we explore the effect of replacing the aggregation weight
$\bm{\alpha}$ with the self-attention aggregation [Shi et al., 2020]
implemented by PyTorch Geometric [Fey and Lenssen, 2019], we name it as “
_SToNe (self-attention)_”.

First of all, we compare the average precision score and AUC score of “
_SToNe_ ” and “ _SToNe (self-attention)_” in Table 6. We observe that “
_SToNe_ ” could achieve superior performance than “ _SToNe (self-attention)_”
on all datasets. In particular, using self-attention with SToNe results in a
larger variance on LastFM and UCI datasets. This is potentially because these
two datasets lack node/edge features and have larger average time-gap (dataset
statistic in Table 3). Since “ _SToNe (self-attention)_” is relying on self-
attention to process the node/edge features, it implicitly assumes node
features exist and are “smooth” at adjacent timestamps. As a result, “ _SToNe
(self-attention)_” could generalize poorly when the assumptions are violated.

Table 6: Comparison on the Average precision and AUC score for link
prediction.

|  | Reddit | Wiki | MOOC | LastFM | GDELT | UCI
---|---|---|---|---|---|---|---
Average Precision |  SToNe (self-attention) | $99.37\pm 0.01$ | $98.99\pm 0.01$ | $99.42\pm 0.01$ | $91.79\pm 4.02$ | $96.20\pm 0.02$ | $90.24\pm 0.73$
SToNe | $\textbf{99.89}\pm 0.00$ | $\textbf{99.85}\pm 0.05$ | $\textbf{99.88}\pm 0.04$ | $\textbf{95.74}\pm 0.13$ | $\textbf{97.40}\pm 0.02$ | $\textbf{94.60}\pm 0.31$
AUC Score |  SToNe (self-attention) | $99.37\pm 0.01$ | $98.92\pm 0.02$ | $99.58\pm 0.01$ | $93.33\pm 4.23$ | $95.56\pm 0.01$ | $88.22\pm 0.98$
SToNe | $\textbf{99.91}\pm 0.00$ | $\textbf{99.85}\pm 0.05$ | $\textbf{99.91}\pm 0.02$ | $\textbf{97.16}\pm 0.06$ | $\textbf{96.96}\pm 0.01$ | $\textbf{93.14}\pm 0.67$

Moreover, to validate the generalization performance, we compare the average
precision score and the generalization gap in Figure 7. The generalization gap
is defined as the absolute difference between the training and validation
average precision scores. Our results show that “ _SToNe_ ” could achieve a
higher average precision score than “ _SToNe (self-attention)_” and has a
smaller generalization gap. This suggests that using aggregation weight
$\bm{\alpha}$ instead of using self-attention in SToNe could lead to better
generalization performance, which is expected because self-attention has
higher model complexity and could be harder to train.

Figure 7: Comparison of the average prevision of the validation set and the
generalization gap of different methods on real-world datasets.

###  B.6 Conduct experiments under the inductive learning setting

Please notice that the discussion in our paper is mainly focusing on
transductive learning. For the completeness, we also report some preliminary
inductive learning results. Our inductive learning setting is the same as [Xu
et al., 2020a]. More specifically, the inductive node sets are all nodes that
does not belonging to the training set nodes, i.e.,
$\mathcal{V}_{\text{inductive}}=\mathcal{V}\setminus\mathcal{V}_{\text{train}}$.
In the inductive setting, negative nodes are only selected from
$\mathcal{V}_{\text{inductive}}$. A test set link is considered as the
inductive link if at least one of its associated nodes belong to
$\mathcal{V}_{\text{inductive}}$, i.e.,
$\mathcal{E}_{\text{inductive}}=\\{(v_{i},v_{j})\leavevmode\nobreak\
|\leavevmode\nobreak\ (v_{i}\in\mathcal{V}_{\text{inductive}}\vee
v_{j}\in\mathcal{V}_{\text{inductive}})\wedge(v_{i},v_{j})\in\mathcal{E}_{\text{test}}\\}$.
Inductive learning performance is evaluated only on the inductive edges.
Dataset statistics are summarized in Table 7.

Table 7: Dataset statistics for inductive learning settings.

| $|\mathcal{V}_{\text{inductive}}|$ | $|\mathcal{E}_{\text{inductive}}|$ | $|\mathcal{E}_{\text{test}}|$
---|---|---|---
Reddit | 134 | 4,704 | 100,867
Wiki | 1,210 | 5,732 | 26,621
MOOC | 342 | 8,645 | 61,763
UCI | 391 | 4,876 | 8,976

We compare with baselines under the inductive learning setting with TGL
framework for a fair comparison. As shown in the table, the performance of our
method also works well on the inductive learning setting. Moreover, by
combining with the results in Table 1, we found that changing from
transductive to inductive learning does not affect our model performance much,
which is potentially because our method has a simpler neural architecture and
a stronger inductive bias for link prediction task.

Table 8: Inductive learning setting average precision.

Inductive average precision | Reddit | Wiki | MOOC | UCI
---|---|---|---|---
SToNe | $99.75\pm 0.00$ | $99.70\pm 0.05$ | $99.81\pm 0.04$ | $90.94\pm 0.33$
TGN | $98.88\pm 0.02$ | $98.88\pm 0.06$ | $99.72\pm 0.08$ | $83.27\pm 2.01$
JODIE | $96.06\pm 0.03$ | $98.59\pm 0.07$ | $86.90\pm 0.06$ | $59.92\pm 3.19$
TGAT | $93.28\pm 0.05$ | $95.82\pm 0.12$ | $95.57\pm 0.11$ | $67.65\pm 0.73$
DySAT | $90.83\pm 0.02$ | $96.45\pm 0.05$ | $98.65\pm 0.19$ | $80.48\pm 2.98$

###  B.7 Preliminary results of node classification

The design of SToNe is mainly focusing on the link prediction. For the
completeness, we also report some preliminary node classification results
using average precision in the table bellow. For this experiment, we use the
identical network structure as the link prediction model. The evaluation
strategy follows the framework in [Zhou et al., 2022]. It is worth noting that
since our paper primarily focused on link prediction task, the performance of
SToNe is not optimal. We plan to make it as an interesting future direction
and enhance model performance by introducing the inductive bias of the node
classification task.

Table 9: Node classification average precision. | SToNe | GraphMixder | JODIE | TGAT | TGN | DySAT
---|---|---|---|---|---|---
Reddit | 70.66 | 70.51 | 70.89 | 63.63 | 64.79 | 62.69
Wiki | 85.78 | 83.93 | 80.63 | 85.30 | 87.13 | 85.30

##  Appendix C Generalization bound of GNN-based method

Recall that we compute the representation of node $v_{i}$ at time $t$ by
applying GNN on the temporal graph that only considers the temporal edges with
timestamp before time $t$. The GNN has trainable weight parameters
$\bm{\theta}=\\{\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(L)}\\}$ and binary hyper-
parameter $\alpha\in\\{0,1\\}$ that controls whether residual connection is
used.

Representation computation. The final prediction on node $v_{i}$ is computed
as $f_{i}(\bm{\theta})=\mathbf{W}^{(L)}\mathbf{h}_{i}^{(L-1)}$, where the
hidden representation $\mathbf{h}_{i}^{(L-1)}$ is computed by

$\displaystyle\mathbf{h}_{i}^{(\ell)}$
$\displaystyle=\sigma\big{(}\mathbf{W}^{(\ell)}\sum\nolimits_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(\ell-1)}\big{)}+\alpha\mathbf{h}_{i}^{(\ell-1)}\in\mathbb{R}^{m},$
$\displaystyle\mathbf{h}_{i}^{(1)}$
$\displaystyle=\sigma\big{(}\mathbf{W}^{(1)}\sum\nolimits_{j\in\mathcal{N}(i)}P_{ij}\mathbf{x}_{j}\big{)}\in\mathbb{R}^{m}.$

Here $\sigma(\cdot)$ is the activation function, $P_{ij}$ is the aggregation
weight used for propagating information from node $v_{j}$ to $v_{i}$, and
$\mathcal{N}(i)$ is the set of all neighbors of node $v_{i}$ in the temporal
graph. For parameter dimension, we have
$\mathbf{W}^{(1)}\in\mathbb{R}^{m\times d}$,
$\mathbf{W}^{(\ell)}\in\mathbb{R}^{m\times m}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ 2\leq\ell\leq L-1$, and
$\mathbf{W}^{(L)}\in\mathbb{R}^{1\times m}$, where $m$ is the hidden dimension
and $d$ is the input dimension.

Gradient computation. The gradient of weight matrix
$\mathbf{W}^{(\ell)},\leavevmode\nobreak\ \forall\ell\in[L-1]$ is computed by

$\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}=\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}\ldots\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}P_{i_{L-1},i_{L-2}}P_{i_{L-2},i_{L-3}}\ldots
P_{i_{\ell+1},i_{\ell}}\cdot\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1}).$

Here $\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1})\in\mathbb{R}^{m\times m}$ is
defined as

$\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1})=\Big{[}\mathbf{W}^{(L)}\Big{(}\mathbf{D}_{i_{L-1}}^{(L-1)}\mathbf{W}^{(L-1)}+\alpha\mathbf{I}_{m}\Big{)}\ldots\Big{(}\mathbf{D}_{i_{\ell+1}}^{(\ell+1)}\mathbf{W}^{(\ell+1)}+\alpha\mathbf{I}_{m}\Big{)}\mathbf{D}_{i_{\ell}}^{(\ell)}\Big{]}^{\top}(\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)})^{\top},$

where
$\mathbf{D}_{i,\ell}^{(\ell)}=\text{diag}(\sigma^{\prime}(\mathbf{z}_{i}^{(\ell)}))\in\mathbb{R}^{m\times
m}$ is a diagonal matrix and $\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}$ is the
aggregation of neighbors’ representations

$\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}=\sum_{i_{\ell-1}\in\mathcal{N}(i_{\ell})}P_{i_{\ell},i_{\ell-1}}\mathbf{h}_{i_{\ell-1}}^{(\ell-1)}\in\mathbb{R}^{m},\leavevmode\nobreak\
\forall\ell\in\\{1,2,\ldots,L-1\\}.$

Finally, the gradient with respect to the final layer weight matrix
$\mathbf{W}^{(L)}$ is computed as

$\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(L)}}=\mathbf{h}^{(L-1)}_{i}.$

Please refer to Appendix G for detailed derivation of the gradients for the
weight parameters in an $L$-layer GNN.

###  C.1 Proof sketch

In the following, we summarize the main steps of proving the generalization
bound of GNN-based methods:

  * •

Firstly, we show in Lemma 1 that if two weight parameters are close to each
other, then the node representation computed on these two parameters are also
close.

  * •

Secondly, we show that if two set of weight parameters are close to each
other, then the neural network outputs computed on these weight parameters are
almost linear in Lemma 2, and the computed loss is almost convex in Lemma 3.

  * •

Thirdly, we show in Lemma 4 that with high probability, the gradient of neural
network can be upper bounded, and this upper bound is dependent on the neural
architecture.

  * •

Thirdly, based on our previous results in Lemma 1, Lemma 3, and Lemma 4, we
show in Lemma 5 that the difference between the cumulative loss over training
to the loss computed on optimal solution could be upper bounded.

  * •

Finally, we show in Lemma 6 that the expected 0-1 error is upper bounded.
Then, by using our definition on the neural tangent random feature and
feature-label alignment, we conclude the proof.

###  C.2 Useful lemmas

For the ease of presentation, we introduce the following two definitions which
will be used when introducing our lemmas.

######  Definition 2 ($\omega$-neighborhood [Cao and Gu, 2019]).

For any
$\widetilde{\bm{\theta}}=\\{\widetilde{\mathbf{W}}^{(1)},\ldots,\widetilde{\mathbf{W}}^{(L)}\\}$,
we define its $\omega$-neighborhood as

$\mathcal{B}(\bm{\theta},\omega)=\\{\widetilde{\bm{\theta}}\leavevmode\nobreak\
|\leavevmode\nobreak\
\|\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)}\|_{\mathrm{F}}\leq\omega,\leavevmode\nobreak\
\ell\in[L]\\}.$

######  Definition 3 (Neural tangent random feature [Cao and Gu, 2019]).

Let $\bm{\theta}_{0}$ be generated via the initialization. We define the
neural tangent random feature function class as

$\mathcal{F}(\bm{\theta}_{0},R)=\\{f(\bm{\theta}_{0})+\langle\nabla_{\theta}f(\bm{\theta}_{0}),\bm{\theta}\rangle\leavevmode\nobreak\
|\leavevmode\nobreak\ \bm{\theta}\in\mathcal{B}(\mathbf{0},Rm^{-1/2})\\},$

where $R>0$ measures the size of the function class and $m$ is the hidden
dimension of the neural network.

In the following, we show that if the input weight parameters
$\bm{\theta}=\\{\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(L)}\\}$ and
$\widetilde{\bm{\theta}}=\\{\widetilde{\mathbf{W}}^{(1)},\ldots,\widetilde{\mathbf{W}}^{(L)}\\}$
are close, the hidden representation of graph neural networks computed on
$\bm{\theta}$ and $\widetilde{\bm{\theta}}$ does not change too much.

######  Lemma 1.

Let $\rho$ be the Lipschitz constant of the activation function and $m$ is the
hidden dimension. Then with $\omega=\mathcal{O}(1/((3\rho+1)\tau)^{(L-1)})$
and assuming $\widetilde{\bm{\theta}}\in\mathcal{B}(\bm{\theta},\omega)$, we
have
$\|\widetilde{\mathbf{h}}_{i}^{(\ell)}-\mathbf{h}_{i}^{(\ell)}\|_{2}=\mathcal{O}(1)$
with probability at least $1-2\ell\exp(-m/2)-\ell\exp(-\Omega(m))$.

Please note that the smaller the distance $\omega$, the closer the
representation
$\|\widetilde{\mathbf{h}}_{i}^{(\ell)}-\mathbf{h}_{i}^{(\ell)}\|_{2}$. In
particular, according to the proof of Lemma 1, we have
$\|\widetilde{\mathbf{h}}_{i}^{(\ell)}-\mathbf{h}_{i}^{(\ell)}\|_{2}\leq\mathcal{O}(\epsilon)$
by selecting $\omega=\mathcal{O}(\epsilon/((3\rho+1)\tau)^{(L-1)})$ for any
small $\epsilon>0$. This conclusion will be later used in Lemma 5.

###### Proof of Lemma 1.

When $\ell=1$, we have for any node $v_{i}\in\mathcal{V}$

$\displaystyle\|\widetilde{\mathbf{h}}_{i}^{(1)}-\mathbf{h}_{i}^{(1)}\|_{2}$
$\displaystyle=\Big{\|}\sigma\Big{(}\widetilde{\mathbf{W}}^{(1)}\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(0)}\Big{)}-\sigma\Big{(}\mathbf{W}^{(1)}\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(0)}\Big{)}\Big{\|}_{2}$
$\displaystyle\underset{(a)}{\leq}\rho\|\widetilde{\mathbf{W}}^{(1)}-\mathbf{W}^{(1)}\|_{2}\Big{\|}\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{x}_{j}\Big{\|}_{2}$
$\displaystyle\underset{(b)}{\leq}\rho\omega\Big{(}\sum_{j\in\mathcal{N}(i)}P_{ij}\|\mathbf{x}_{j}\|_{2}\Big{)}$
$\displaystyle\underset{(c)}{\leq}\rho\tau\cdot\omega=\mathcal{O}(1),$

where inequality (a) is due to the Lipschitz continuity of activation
function, inequality (b) is due to the $\omega$-neighborhood definition
$\widetilde{\mathbf{W}}\in\mathcal{B}(\mathbf{W},\omega)$, and inequality (c)
is due to Assumption 1 and Assumption 3.

Similarly, when $\ell\in\\{2,\ldots,L-1\\}$, we have $\forall i\in\mathcal{V}$

$\displaystyle\|\widetilde{\mathbf{h}}_{i}^{(\ell)}-\mathbf{h}_{i}^{(\ell)}\|_{2}$
$\displaystyle\underset{(a)}{\leq}\Big{\|}\sigma\Big{(}\widetilde{\mathbf{W}}^{(\ell)}\sum_{j\in\mathcal{N}(i)}P_{ij}\widetilde{\mathbf{h}}_{j}^{(\ell-1)}\Big{)}-\sigma\Big{(}\mathbf{W}^{(\ell)}\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(\ell-1)}\Big{)}\Big{\|}_{2}+\alpha\|\widetilde{\mathbf{h}}_{i}^{(\ell-1)}-\mathbf{h}_{i}^{(\ell-1)}\|_{2}$
$\displaystyle\underset{(b)}{\leq}\rho\|\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)}\|_{2}\Big{\|}\sum_{j\in\mathcal{N}(i)}P_{ij}\widetilde{\mathbf{h}}_{j}^{(\ell-1)}\Big{\|}_{2}+\rho\|\mathbf{W}^{(\ell)}\|_{2}\Big{(}\sum_{j\in\mathcal{N}(i)}P_{ij}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2}\Big{)}$
$\displaystyle\quad+\|\widetilde{\mathbf{h}}_{i}^{(\ell-1)}-\mathbf{h}_{i}^{(\ell-1)}\|_{2}$
$\displaystyle\underset{(c)}{\leq}\rho\omega\Big{(}\Big{\|}\sum_{j\in\mathcal{N}(i)}P_{ij}\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(\ell-1)}\Big{\|}_{2}+\Big{\|}\sum_{j\in\mathcal{N}(i)}P_{ij}\mathbf{h}_{j}^{(\ell-1)}\Big{\|}_{2}\Big{)}$
$\displaystyle\quad+\rho\|\mathbf{W}^{(\ell)}\|_{2}\Big{(}\sum_{j\in\mathcal{N}(i)}P_{ij}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2}\Big{)}+\|\widetilde{\mathbf{h}}_{i}^{(\ell-1)}-\mathbf{h}_{i}^{(\ell-1)}\|_{2}$
$\displaystyle\leq\rho\omega\tau\cdot\max_{j\in\\{i\\}\cup\mathcal{N}(i)}\|\mathbf{h}_{j}^{(\ell-1)}\|_{2}+\Big{(}\rho\tau(\|\mathbf{W}^{(\ell)}\|_{2}+\omega)+1\Big{)}\cdot\max_{j\in\\{i\\}\cup\mathcal{N}(i)}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2},$

where the inequality (a) and (c) are due to
$\|\mathbf{A}+\mathbf{B}\|_{2}\leq\|\mathbf{A}\|_{2}+\|\mathbf{B}\|_{2}$, the
inequality (b) is due to the Lipschitz continuity of activation function and
$\widetilde{\mathbf{W}}\in\mathcal{B}(\mathbf{W},\omega)$.

By Proposition 19, we know that with probability at least $1-2\exp(-m/2)$ we
have $\|\mathbf{W}^{(\ell)}\|_{2}\leq 3$ for all $\ell\in[L-1]$.

By Lemma 21, we know that with probability at least $1-\exp(-\Omega(m))$ we
have $\|\mathbf{h}_{i}^{(\ell)}\|_{2}=\Theta(1)$.

Then, we have with probability at least
$1-2\ell\exp(-m/2)-\ell\exp(-\Omega(m))$ for any $i\in\mathcal{V}$

$\displaystyle\|\widetilde{\mathbf{h}}_{i}^{(\ell)}-\mathbf{h}_{i}^{(\ell)}\|_{2}$
$\displaystyle\leq\Big{(}\rho\tau(3+\omega)+1\Big{)}\cdot\max_{j\in\\{i\\}\cup\mathcal{N}(i)}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2}+\rho\omega\tau\cdot\max_{j\in\\{i\\}\cup\mathcal{N}(i)}\|\mathbf{h}_{j}^{(\ell-1)}\|_{2}$
$\displaystyle\leq\Big{(}\rho\tau(3+\omega)+1\Big{)}\cdot\max_{j\in\mathcal{V}}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2}+\rho\omega\tau\cdot\max_{j\in\mathcal{V}}\|\mathbf{h}_{j}^{(\ell-1)}\|_{2}$
$\displaystyle\underset{(a)}{=}\Big{(}\rho\tau(3+\omega)+1\Big{)}\cdot\max_{j\in\mathcal{V}}\|\widetilde{\mathbf{h}}_{j}^{(\ell-1)}-\mathbf{h}_{j}^{(\ell-1)}\|_{2}+\rho\omega\tau\cdot\Theta(1)$
$\displaystyle\leq\rho\omega\tau\cdot\Theta(1)\cdot\frac{((3+\omega)\rho\tau+1)^{\ell-1}-1}{(3+\omega)\rho\tau}+\rho\omega\tau\cdot((3+\omega)\rho\tau+1)^{\ell-1}$
$\displaystyle\leq((3+\omega)\rho\tau+1)^{\ell-1}\Big{(}\frac{\rho\omega\tau\cdot\Theta(1)}{(3+\omega)\rho\tau}+\rho\omega\tau\Big{)}$
$\displaystyle\leq((3+\omega)\rho\tau+1)^{\ell-1}\cdot\omega\cdot\Big{(}\Theta(1)+\rho\tau\Big{)},$

where the equality (a) is due to
$\max_{j\in\mathcal{V}}\|\mathbf{h}_{j}^{(\ell)}\|_{2}=\Theta(1)$.

By setting $\omega=1/((3\rho+1)\tau)^{L-1}$ we have the above equation upper
bounded by $\mathcal{O}(1)$. ∎

Then, in the next lemma, we show that if the initialization of two set of
weight parameters $\bm{\theta}=\\{\mathbf{W}^{(1)},\ldots,\mathbf{W}^{(L)}\\}$
and
$\widetilde{\bm{\theta}}=\\{\widetilde{\mathbf{W}}^{(1)},\ldots,\widetilde{\mathbf{W}}^{(L)}\\}$
are close, the neural network output $f_{i}(\bm{\theta})$ will be almost
linear with respect to its weight parameters.

######  Lemma 2.

Let
$\bm{\theta},\widetilde{\bm{\theta}}\in\mathcal{B}(\bm{\theta}_{0},\omega)$
with $\omega=\mathcal{O}\Big{(}1/((3\rho+1)\tau)^{(L-1)}\Big{)}$. Then, for
any node $v_{i}\in\mathcal{V}$ in the graph, with probability at least
$1-2(L-1)\exp(-m/2)-L\exp(-\Omega(m))-2/m$, we have
$\epsilon_{\text{lin}}:=|f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})-\langle\nabla
f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\rangle|=\mathcal{O}(1),$
where $f_{i}(\bm{\theta})$ is the prediction on the $L$-hop subgraph centered
on the root node $v_{i}$.

Please note that the smaller the distance $\omega$, the more the model output
close to linear. In particular, according to the proof of Lemma 2 and the
proof of Lemma 1, we have $\epsilon_{\text{lin}}\leq\mathcal{O}(\epsilon)$ by
selecting $\omega=\mathcal{O}(\epsilon/((3\rho+1)\tau)^{(L-1)})$ for any small
$\epsilon>0$. This conclusion will be later used in Lemma 5.

###### Proof of Lemma 2.

According to the forward and backward propagation rules as we recapped at the
beginning of this section, we have

$\displaystyle\left|f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})-\langle\nabla
f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\rangle\right|$
$\displaystyle=\left|f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})-\sum_{\ell=1}^{L}\Big{\langle}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}},\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)}\Big{\rangle}\right|$
$\displaystyle\leq\left|f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})-\Big{\langle}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(L)}},\widetilde{\mathbf{W}}^{(L)}-\mathbf{W}^{(L)}\Big{\rangle}\right|+\sum_{\ell=1}^{L-1}\left|\Big{\langle}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}},\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)}\Big{\rangle}\right|$
$\displaystyle=\left|\widetilde{\mathbf{W}}^{(L)}(\widetilde{\mathbf{h}}_{i}^{(L-1)}-\mathbf{h}_{i}^{(L-1)})^{\top}\right|$
$\displaystyle+\sum_{\ell=1}^{L-1}\Big{[}\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}\ldots\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}P_{i_{L-1},i_{L-2}}P_{i_{L-2},i_{L-3}}\ldots
P_{i_{\ell+1},i_{\ell}}$
$\displaystyle\quad\Big{|}\mathbf{W}^{(L)}\Big{(}\mathbf{D}_{i_{L-1}}^{(L-1)}\mathbf{W}^{(L-1)}+\alpha_{L-2}\mathbf{I}_{m}\Big{)}\ldots\Big{(}\mathbf{D}_{i_{\ell+1}}^{(\ell+1)}\mathbf{W}^{(\ell+1)}+\alpha_{\ell}\mathbf{I}_{m}\Big{)}\mathbf{D}_{i_{\ell}}^{(\ell)}(\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)})\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}\Big{|}\Big{]}$
$\displaystyle\leq\|\widetilde{\mathbf{W}}^{(L)}\|_{2}\Big{\|}\widetilde{\mathbf{h}}_{i}^{(L-1)}-\mathbf{h}_{i}^{(L-1)}\Big{\|}_{2}$
$\displaystyle+\sum_{\ell=1}^{L-1}\Big{[}\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}\ldots\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}P_{i_{L-1},i_{L-2}}P_{i_{L-2},i_{L-3}}\ldots
P_{i_{\ell+1},i_{\ell}}$
$\displaystyle\qquad\underbrace{\|\mathbf{W}^{(L)}\|_{2}\underbrace{\Big{\|}\Big{(}\mathbf{D}_{i_{L-1}}^{(L-1)}\mathbf{W}^{(L-1)}+\alpha_{L-2}\mathbf{I}_{m}\Big{)}\ldots\Big{(}\mathbf{D}_{i_{\ell+1}}^{(\ell+1)}\mathbf{W}^{(\ell+1)}+\alpha_{\ell}\mathbf{I}_{m}\Big{)}\mathbf{D}_{i_{\ell}}^{(\ell)}\Big{\|}_{2}}_{(a)}\|\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}\|_{2}\Big{\|}\widetilde{\mathbf{W}}^{(\ell)}-\mathbf{W}^{(\ell)}\Big{\|}_{2}}_{(b)}\Big{]},$

where
$\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}=\sum_{i_{\ell-1}\in\mathcal{N}(i_{\ell})}P_{{i_{\ell},i_{\ell-1}}}\mathbf{h}_{i_{\ell-1}}^{(\ell-1)}$
and it has bounded $\ell_{2}$-norm

$\|\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}\|_{2}\leq\tau\max_{i_{\ell-1}\in\mathcal{N}(i_{\ell})}\|\mathbf{h}_{i_{\ell-1}}^{(\ell-1)}\|_{2}=\tau\Theta(1).$

Since the derivative of activation function is bounded, we have
$\|\mathbf{D}_{i,\ell}\|_{2}\leq\rho$ and
$\|\mathbf{D}_{i}^{(\ell)}\mathbf{W}^{(\ell)}+\alpha_{\ell}\mathbf{I}_{m}\|_{2}\leq
3\rho+1$.

Therefore, we know that the term (a) in the above equation could be upper
bounded by $\rho(3\rho+1)^{L-1-\ell}$.

Besides, we know that $\|\mathbf{W}^{(L)}\|_{2}\leq\sqrt{2}$ according to
Lemma 20.

Finally, by plugging the results back, we can upper bound the term (b) in the
above equation by

$(b)\leq C_{\ell}=\sqrt{2}\rho(3\rho+1)^{L-\ell-1}\tau\Theta(1)\omega$

and

$\displaystyle\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}\ldots\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}\Big{(}P_{i_{L-1},i_{L-2}}P_{i_{L-2},i_{L-3}}\ldots
P_{i_{\ell+1},i_{\ell}}C_{\ell}\Big{)}$
$\displaystyle=\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}P_{i_{L-1},i_{L-2}}\Big{(}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}P_{i_{L-2},i_{L-3}}\ldots\Big{(}\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}P_{i_{\ell+1},i_{\ell}}C_{\ell}\Big{)}\ldots\Big{)}$
$\displaystyle=\sqrt{2}\rho(3\rho+1)^{L-\ell-1}\tau^{L-\ell}\Theta(1)\omega.$

As a result, we have

$\displaystyle|f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})-\langle\nabla
f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\rangle|$
$\displaystyle\leq(\sqrt{2}+\omega)\cdot\mathcal{O}(1)+\sqrt{2}\rho\omega\cdot\Theta(1)\cdot\sum_{\ell=1}^{L-1}(3\rho+1)^{L-\ell-1}\tau^{L-\ell}$
$\displaystyle\leq(\sqrt{2}+\omega)\cdot\mathcal{O}(1)+\sqrt{2}\rho\omega\cdot\Theta(1)\cdot\Big{(}\frac{((3\rho+1)\tau)^{L-1}-1}{3\rho}\Big{)}$
$\displaystyle\leq\mathcal{O}(1),$

where the last inequality holds by selecting $\omega=1/((3\rho+1)\tau)^{L-1}$.
∎

Let us define the logistic loss as

$\text{loss}_{i}(\bm{\theta})=\psi(y_{i}f_{i}(\bm{\theta})),\psi(x)=\log(1+\exp(-x)).$

Then, the following lemma shows that $\text{loss}_{i}(\bm{\theta})$ is almost
a convex function of $\bm{\theta}$ for any $v_{i}\in\mathcal{V}$ if the
initialization of two set of parameters are close to each other.

######  Lemma 3.

Let
$\bm{\theta},\widetilde{\bm{\theta}}\in\mathcal{B}(\bm{\theta}_{0},\omega)$
with $\omega=1/((3\rho+1)\tau)^{L-1}$ for any $v_{i}\in\mathcal{V}$, it holds
that
$\text{loss}_{i}(\widetilde{\bm{\theta}})\geq\text{loss}_{i}(\bm{\theta})+\langle\nabla_{\theta}\text{loss}_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\rangle-\epsilon_{\text{lin}}$
with probability at least $1-2(L-1)\exp(-m/2)-L\exp(-\Omega(m))-2/m$, where
$\epsilon_{\text{lin}}=\Big{|}\Big{\langle}\nabla_{\theta}f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}-f_{i}(\widetilde{\bm{\theta}})+f_{i}(\bm{\theta})\Big{|}=\mathcal{O}(1)$
according to Lemma 2.

###### Proof of Lemma 3.

The proof follows the proof of Lemma 9 in [Zhu et al., 2022].

By the convexity of $\psi(\cdot)$, we know that
$\psi(b)\geq\psi(a)+\psi^{\prime}(a)(b-a)$. Therefore, we have

$\displaystyle\text{loss}_{i}(\widetilde{\bm{\theta}})-\text{loss}_{i}(\bm{\theta})$
$\displaystyle=\psi(y_{i}f_{i}(\widetilde{\bm{\theta}}))-\psi(y_{i}f_{i}(\bm{\theta}))$
$\displaystyle\geq\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))\Big{(}y_{i}f_{i}(\widetilde{\bm{\theta}})-y_{i}f_{i}(\bm{\theta})\Big{)}$
$\displaystyle=\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))y_{i}\Big{(}f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})\Big{)}.$

By using the chain rule, we have

$\Big{\langle}\nabla_{\theta}\text{loss}_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}=\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))\cdot
y_{i}\Big{\langle}\nabla_{\theta}f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}.$

By combining the above equations, we have

$\displaystyle\text{loss}_{i}(\widetilde{\bm{\theta}})-\text{loss}_{i}(\bm{\theta})$
$\displaystyle\geq\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))y_{i}\Big{(}f_{i}(\widetilde{\bm{\theta}})-f_{i}(\bm{\theta})\Big{)}$
$\displaystyle=\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))y_{i}\Big{\langle}\nabla_{\theta}f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}-\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))y_{i}\Big{(}\Big{\langle}\nabla_{\theta}f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}-f_{i}(\widetilde{\bm{\theta}})+f_{i}(\bm{\theta})\Big{)}$
$\displaystyle\geq\Big{\langle}\nabla_{\theta}\text{loss}_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}-\Big{|}\Big{\langle}\nabla_{\theta}f_{i}(\bm{\theta}),\widetilde{\bm{\theta}}-\bm{\theta}\Big{\rangle}-f_{i}(\widetilde{\bm{\theta}})+f_{i}(\bm{\theta})\Big{|},$

where the inequality is due to
$|\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))y_{i}|\leq 1$. ∎

Moreover, by the gradient computation, we know that the gradient of the neural
network function can be upper bounded.

######  Lemma 4.

For any $v_{i}\in\mathcal{V}$ with probability at least
$1-2(L-\ell)\exp(-m/2)-\ell\exp(-\Omega(m))-2/m$, it holds that
$\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2},\Big{\|}\frac{\partial\text{loss}_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}\leq\Theta\Big{(}((3\rho+1)\tau)^{L-\ell}\Big{)}$

###### Proof of Lemma 4.

For $\ell=L$, we have for any $v_{i}\in\mathcal{V}$

$\displaystyle\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(L)}}\Big{\|}_{2}$
$\displaystyle=\|\mathbf{h}_{i}^{(L-1)}\|_{2}$
$\displaystyle\leq\max_{j\in\mathcal{V}}\|\mathbf{h}_{j}^{(L-1)}\|_{2}=\Theta(1).$

For $\ell\in[L-1]$, we have for any $v_{i}\in\mathcal{V}$

$\displaystyle\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}\leq\sum_{i_{L-2}\in\mathcal{N}(i_{L-1})}\sum_{i_{L-3}\in\mathcal{N}(i_{L-2})}\ldots\sum_{i_{\ell}\in\mathcal{N}(i_{\ell+1})}P_{i_{L-1},i_{L-2}}P_{i_{L-2},i_{L-3}}\ldots
P_{i_{\ell+1},i_{\ell}}\Big{\|}\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1})\Big{\|}_{2}.$

To upper bound
$\Big{\|}\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1})\Big{\|}_{2}$, we have

$\displaystyle\Big{\|}\mathbf{G}^{\ell}(i_{\ell},\ldots,i_{L-1})\Big{\|}_{2}$
$\displaystyle=\Big{\|}\Big{[}\mathbf{W}^{(L)}\Big{(}\mathbf{D}_{i_{L-1}}^{(L-1)}\mathbf{W}^{(L-1)}+\alpha_{L-2}\mathbf{I}_{m}\Big{)}\ldots\Big{(}\mathbf{D}_{i_{\ell+1}}^{(\ell+1)}\mathbf{W}^{(\ell+1)}+\alpha_{\ell}\mathbf{I}_{m}\Big{)}\mathbf{D}_{i_{\ell}}^{(\ell)}\Big{]}^{\top}(\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)})^{\top}\Big{\|}_{2}$
$\displaystyle\leq\Big{\|}\mathbf{W}^{(L)}\Big{\|}_{2}\Big{\|}\Big{(}\mathbf{D}_{i_{L-1}}^{(L-1)}\mathbf{W}^{(L-1)}+\alpha_{L-2}\mathbf{I}_{m}\Big{)}\ldots\Big{(}\mathbf{D}_{i_{\ell+1}}^{(\ell+1)}\mathbf{W}^{(\ell+1)}+\alpha_{\ell}\mathbf{I}_{m}\Big{)}\mathbf{D}_{i_{\ell}}^{(\ell)}\Big{\|}_{2}\Big{\|}\tilde{\mathbf{z}}_{i_{\ell}}^{(\ell-1)}\Big{\|}_{2}$
$\displaystyle\leq\sqrt{2}\rho\tau(3\rho+1)^{L-\ell-1}\Theta(1),$

where the last inequality follows the proof of Lemma 2.

By combining the above results, we have

$\displaystyle\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}$
$\displaystyle\leq\sqrt{2}\rho((3\rho+1)\tau)^{L-\ell-1}\Theta(1)$
$\displaystyle\leq\Theta((3\rho+1)\tau)^{L-\ell}.$

Moreover, the above inequality also implies

$\displaystyle\Big{\|}\frac{\partial\text{loss}_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}$
$\displaystyle=|\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))\cdot
y_{i}|\cdot\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}$
$\displaystyle\leq\Big{\|}\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2},$

where the last inequality is due to
$|\psi^{\prime}(y_{i}f_{i}(\bm{\theta}))\cdot y_{i}|\leq 1$. ∎

In the following, we show that the cumulative loss can be upper bounded under
small changes on the weight parameters.

######  Lemma 5.

For any $\epsilon,\delta,R>0$, there exists
$m^{\star}=\mathcal{O}\left(\frac{((3\rho+1)\tau)^{4(L-1)}L^{2}R^{4}}{4\epsilon^{4}}\right)\log(1/\delta),$
such that if $m\geq m^{\star}$, then with probability at least $1-\delta$ over
the randomness of $\bm{\theta}_{0}$, for any
$\bm{\theta}_{\star}\in\mathcal{B}(\bm{\theta}_{0},Rm^{-1/2})$, with
$\eta=\frac{\epsilon}{mL((3\rho+1)\tau)^{2(L-1)}}$ and
$N=\frac{L^{2}R^{2}((3\rho+1)\tau)^{2(L-1)}}{2\epsilon^{2}}$, the cumulative
loss can be upper bounded by
$\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{i-1})\leq\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}^{\star})+\mathcal{O}(\epsilon).$

###### Proof of Lemma 5.

Let us define $\bm{\theta}^{\star}$ as the optimal solution that could
minimize the cumulative loss over $N$ epochs, where at each epoch only a
single data point is used as defined in Algorithm 1

$\bm{\theta}^{\star}=\underset{\bm{\theta}\in\mathcal{B}(\bm{\theta}_{0},\omega)}{\arg\min}\sum\nolimits_{i=1}^{N}\text{loss}_{i}(\bm{\theta}).$

Without loss of generality, let us assume the epoch loss
$\text{loss}_{i}(\bm{\theta}_{i-1})$ is computed on the $i$-th node. Then, in
the following, we try to show
$\bm{\theta}_{0},\bm{\theta}_{1},\ldots,\bm{\theta}_{N-1}\in\mathcal{B}(\bm{\theta}_{0},\omega)$,
where $\omega=\epsilon/((3\rho+1)\tau)^{L-1}$

First of all, it is clear that
$\bm{\theta}_{0}\in\mathcal{B}(\bm{\theta}_{0},\omega)$. Then, to show
$\bm{\theta}_{n}\in\mathcal{B}(\bm{\theta}_{0},\omega)$ for any
$n\in\\{1,\ldots,N-1\\}$, we use our previous conclusion on the upper bound of
gradient
$\|\partial\text{loss}_{i}(\bm{\theta})/\partial\mathbf{W}^{(\ell)}\|\leq\Theta((3\rho+1)\tau)^{L-\ell}$
in Lemma 4 and have

$\displaystyle\|\mathbf{W}^{(\ell)}_{n}-\mathbf{W}^{(\ell)}_{0}\|_{2}$
$\displaystyle\leq\sum_{i=1}^{n}\|\mathbf{W}^{(\ell)}_{i}-\mathbf{W}^{(\ell)}_{i-1}\|_{2}$
$\displaystyle=\sum_{i=1}^{n}\eta\Big{\|}\frac{\partial\text{loss}_{i-1}(\bm{\theta}_{i-1})}{\partial\mathbf{W}^{(\ell)}_{i-1}}\Big{\|}_{2}$
$\displaystyle\leq\Theta\Big{(}\eta N((3\rho+1)\tau)^{L-\ell}\Big{)}.$

By plugging in the choice of
$\eta=\frac{\epsilon}{mL((3\rho+1)\tau)^{2(L-1)}}$ and
$N=\frac{L^{2}R^{2}((3\rho+1)\tau)^{2(L-1)}}{2\epsilon^{2}}$, we have

$\displaystyle\|\mathbf{W}^{(\ell)}_{n}-\mathbf{W}^{(\ell)}_{0}\|_{2}$
$\displaystyle\leq\Theta\Big{(}\eta N((3\rho+1)\tau)^{L-\ell-1}\Big{)}$
$\displaystyle=\Theta\Big{(}\frac{\epsilon}{mL((3\rho+1)\tau)^{2(L-1)}}\frac{L^{2}R^{2}((3\rho+1)\tau)^{2(L-1)}}{2\epsilon^{2}}((3\rho+1)\tau)^{L-\ell}\Big{)}$
$\displaystyle=\Theta\Big{(}\frac{LR^{2}}{2m\epsilon}((3\rho+1)\tau)^{L-\ell}\Big{)}.$

After changing norm from $\ell_{2}$-norm to Frobenius-norm, we have

$\|\mathbf{W}^{(\ell)}_{n}-\mathbf{W}^{(\ell)}_{0}\|_{\mathrm{F}}\leq\sqrt{m}\times\Theta\Big{(}\frac{LR^{2}}{2m\epsilon}((3\rho+1)\tau)^{L-\ell}\Big{)},$
(2)

By plugging in the selection of hidden dimension $m$, we have

$\|\mathbf{W}^{(\ell)}_{n}-\mathbf{W}^{(\ell)}_{0}\|_{\mathrm{F}}\leq\frac{\epsilon}{((3\rho+1)\tau)^{L-1}},$

which means
$\bm{\theta}_{0},\bm{\theta}_{1},\ldots,\bm{\theta}_{N-1}\in\mathcal{B}(\bm{\theta}_{0},\omega)$
for $\omega=\epsilon/((3\rho+1)\tau)^{L-1}$

Then, our next step is to bound
$\text{loss}_{i}(\bm{\theta}_{i})-\text{loss}_{i}(\bm{\theta}_{\star})$. By
Lemma 3, we know that

$\displaystyle\text{loss}_{i+1}(\bm{\theta}_{i})-\text{loss}_{i+1}(\bm{\theta}_{\star})$
$\displaystyle\leq\Big{\langle}\nabla_{\theta}\text{loss}_{i+1}(\bm{\theta}_{i}),\bm{\theta}_{i}-\bm{\theta}_{\star}\Big{\rangle}+\epsilon_{\text{lin}}$
$\displaystyle=\sum_{\ell=1}^{L}\Big{\langle}\frac{\partial\text{loss}_{i+1}(\bm{\theta}_{i})}{\partial\mathbf{W}^{(\ell)}},\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\Big{\rangle}+\epsilon_{\text{lin}}$
$\displaystyle=\frac{1}{\eta}\sum_{\ell=1}^{L}\Big{\langle}\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{i+1}^{(\ell)},\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\Big{\rangle}+\epsilon_{\text{lin}}$
$\displaystyle\leq\frac{1}{2\eta}\sum_{\ell=1}^{L}\Big{(}\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{i+1}^{(\ell)}\|_{\mathrm{F}}^{2}+\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}-\|\mathbf{W}_{i+1}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}\Big{)}+\epsilon_{\text{lin}}.$

Then, our next step is to upper bound each term on the right hand size of
inequality.

(1) According to the proof of Lemma 2, we have
$\epsilon_{\text{lin}}\leq\mathcal{O}(\epsilon)$ by selecting
$\omega=\mathcal{O}(\epsilon/((3\rho+1)\tau)^{(L-1)})$.

(2) Recall that
$\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{i+1}^{(\ell)}\|_{\mathrm{F}}^{2}$ could
be upper bounded by

$\displaystyle\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{i+1}^{(\ell)}\|_{\mathrm{F}}^{2}$
$\displaystyle\leq\eta^{2}\Big{\|}\frac{\partial\text{loss}_{i+1}(\bm{\theta}_{i})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{\mathrm{F}}^{2}$
$\displaystyle\leq\eta^{2}m\Big{\|}\frac{\partial\text{loss}_{i+1}(\bm{\theta}_{i})}{\partial\mathbf{W}^{(\ell)}}\Big{\|}_{2}^{2}$
$\displaystyle\leq\Theta\Big{(}m\eta^{2}((3\rho+1)\tau)^{2(L-\ell)}\Big{)}.$

(3) By finite sum
$\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}-\|\mathbf{W}_{i+1}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}$
from $i=0$ to $N-1$, we have

$\displaystyle\frac{1}{N}\sum_{i=0}^{N-1}(\|\mathbf{W}_{i}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}-\|\mathbf{W}_{i+1}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2})$
$\displaystyle=\frac{1}{N}\|\mathbf{W}_{0}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}\underbrace{-\frac{1}{N}\|\mathbf{W}_{N+1}^{(\ell)}-\mathbf{W}_{\star}^{(\ell)}\|_{\mathrm{F}}^{2}}_{\leq
0}$ $\displaystyle\leq\frac{R^{2}}{mN},$

where the inequality is due to
$\bm{\theta}_{\star}\in\mathcal{B}(\bm{\theta}_{0},Rm^{-1/2})$.

Finally, by combining the results above, we have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{i-1})-\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{\star})$
$\displaystyle\leq\frac{R^{2}L}{2m\eta
N}+\sum_{\ell=1}^{L}\Theta\Big{(}\frac{m\eta}{2}(3\rho\tau+1)^{2(L-1)}\Big{)}+\mathcal{O}(\epsilon).$

By selecting $\eta=\frac{\epsilon}{mL((3\rho+1)\tau)^{2(L-1)}}$ and
$N=\frac{L^{2}R^{2}((3\rho+1)\tau)^{2(L-1)}}{2\epsilon^{2}}$, we have

$\displaystyle\frac{R^{2}L}{2m\eta N}$
$\displaystyle=\frac{R^{2}L}{2m}\times\frac{mL((3\rho+1)\tau)^{2(L-1)}}{\epsilon}\times\frac{2\epsilon^{2}}{L^{2}R^{2}(3\rho\tau+1)^{2(L-1)}}=\epsilon,$
$\displaystyle\sum_{\ell=1}^{L}\Theta\Big{(}\frac{m\eta}{2}((3\rho+1)\tau)^{2(L-\ell)}\Big{)}$
$\displaystyle\leq
L\times\Theta\Big{(}m\eta\cdot((3\rho+1)\tau)^{2(L-1)}\Big{)}$
$\displaystyle=\Theta\Big{(}m\cdot\frac{\epsilon}{mL((3\rho+1)\tau)^{2(L-1)}}((3\rho+1)\tau)^{2(L-1)}\Big{)}=\Theta(\epsilon).$

Therefore, combining the results above, we have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{i-1})-\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{\star})$
$\displaystyle\leq\mathcal{O}(\epsilon)$

∎

By plugging in the selection of
$\epsilon=\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{2N}}$ to Lemma 5, we have

######  Corollary 1.

For any $\delta>0$ and $R>0$, there exists
$m^{\star}=\mathcal{O}\left(N^{2}/L^{2}\right)\log(1/\delta),$ such that if
$m\geq m^{\star}$, then with probability at least $1-\delta$ over the
randomness of $\bm{\theta}_{0}$, for any
$\bm{\theta}_{\star}\in\mathcal{B}(\bm{\theta}_{0},Rm^{-1/2})$, with the
selection of learning rate $\eta=\frac{R}{m\sqrt{2N}((3\rho+1)\tau)^{(L-1)}}$,
the cumulative loss can be upper bounded by
$\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}_{i-1})\leq\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}^{\star})+\mathcal{O}\left(\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\right).$

In the following, we present the expected 0-1 error bound of multi-layer GNN,
which consists of two terms: (1) the expected 0-1 error with the neural
tangent random feature function and (2) the standard large-deviation error
term.

######  Lemma 6.

For any $\delta\in(0,1/e]$ and $R>0$, there exists
$m^{\star}=\mathcal{O}\left(N^{2}/L^{2}\right)\log(1/\delta),$ such that if
$m\geq m^{\star}$, then with probability at least $1-\delta$ over the
randomness of $\bm{\theta}_{0}$, with the selection of learning rate
$\eta=\frac{R}{m\sqrt{2N}((3\rho+1)\tau)^{(L-1)}}$, we have
$\displaystyle\mathbb{E}\left[\text{loss}_{N}^{0-1}(\widetilde{\bm{\theta}})|\mathcal{D}_{1}^{N-1}\right]$
$\displaystyle\leq\frac{4}{N}\inf_{f\in\mathcal{F}(\bm{\theta}_{0},R)}\sum_{i=1}^{N}\psi(y_{i}f_{i})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}+\mathcal{O}\Big{(}\sqrt{\frac{\log(1/\delta)}{N}}\Big{)}+\Delta,$
where $\mathcal{F}(\bm{\theta}_{0},R)$ is the neural tangent random feature
function class, $\widetilde{\bm{\theta}}$ is uniformly selected from
$\\{\bm{\theta}_{0},\ldots,\bm{\theta}_{N-1}\\}$, $\mathcal{D}_{1}^{N-1}$ is
the sequence of data points sampled before the $N$-th iteration, and the
expectation is computed on the uniform selection of weight parameters
$\widetilde{\bm{\theta}}$ and the condition sampling of $N$-th iteration data
examples.

###### Proof of Lemma 6.

The proof is based on the proof of Theorem 3.3 in [Cao and Gu, 2019].

Let us recall that the 0-1 loss is defined as

$\text{loss}^{0-1}_{i}(\bm{\theta})=\mathbf{1}\\{y_{i}f_{i}(\bm{\theta})<0\\},\leavevmode\nobreak\
\forall i\in[N]$

Since the cross entropy loss $\psi(\cdot)$ satisfies $\mathbf{1}\\{z\leq
0\\}\leq 4\psi(z)$, we have $\text{loss}^{0-1}_{i}(\bm{\theta})\leq
4\text{loss}_{i}(\bm{\theta})$. Then, we have with probability at least
$1-\delta$

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}^{0-1}(\bm{\theta}_{i-1})$
$\displaystyle\leq\frac{4}{N}\sum_{i=1}^{N}\text{loss}_{i}(\bm{\theta}^{\star})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}$
$\displaystyle\leq\frac{4}{N}\sum_{i=1}^{N}\psi(y_{i}f_{i}(\bm{\theta}_{\star}))+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}$

Let us define
$F_{i}(\bm{\theta}_{0},\bm{\theta}^{\star}):=f_{i}(\bm{\theta}_{0})+\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}^{\star}-\bm{\theta}_{0}\rangle$.
Since $\psi(\cdot)$ is $1$-Lipschitz continuous, we have

$\displaystyle\psi(y_{i}f_{i}(\bm{\theta}))-\psi(y_{i}F_{i}(\bm{\theta}_{0},\bm{\theta}^{\star}))$
$\displaystyle\leq
y_{i}\Big{(}f_{i}(\bm{\theta})-F_{i}(\bm{\theta}_{0},\bm{\theta}^{\star})\Big{)}$
$\displaystyle=y_{i}\Big{(}f_{i}(\bm{\theta})-f_{i}(\bm{\theta}_{0})+\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}_{0}-\bm{\theta}^{\star}\rangle\Big{)}\leq\mathcal{O}(1),$

where the last inequality is due to Lemma 2.

By combining the results above, we have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\text{loss}_{i}^{0-1}(\bm{\theta}_{i-1})$
$\displaystyle\leq\frac{4}{N}\sum_{i=1}^{N}\psi(y_{i}F_{i}(\bm{\theta}_{0},\bm{\theta}^{\star}))+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}$
$\displaystyle\leq\frac{4}{N}\inf_{f\in\mathcal{F}(\bm{\theta}_{0},R)}\sum_{i=1}^{N}\psi(y_{i}f_{i})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)},$

where the second inequality is due to the definition of neural tangent random
feature

$\mathcal{F}(\bm{\theta}_{0},R)=\\{f(\bm{\theta}_{0})+\langle\nabla_{\theta}f(\bm{\theta}_{0}),\bm{\theta}\rangle\leavevmode\nobreak\
|\leavevmode\nobreak\ \bm{\theta}\in\mathcal{B}(\mathbf{0},Rm^{-1/2})\\}$

and $\bm{\theta}^{\star}\in\mathcal{B}(\bm{\theta}_{0},Rm^{-1/2})$.

By Proposition 1, we have with probability at least $1-\delta$,

$\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\left[\text{loss}_{i}^{0-1}(\bm{\theta}_{i-1})|\mathcal{D}_{1}^{i-1}\right]\leq\frac{4}{N}\inf_{f\in\mathcal{F}(\bm{\theta}_{0},R)}\sum_{i=1}^{N}\psi(y_{i}f_{i})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}+\mathcal{O}\Big{(}\sqrt{\frac{\log(1/\delta)}{N}}\Big{)}$

By using Assumption 4, we have

$\displaystyle\mathbb{E}\left[\text{loss}_{N}^{0-1}(\bm{\theta})|\mathcal{D}_{1}^{N-1}\right]$
$\displaystyle=\frac{1}{N}\sum_{i=1}^{N}\mathbb{E}\left[\text{loss}_{N}^{0-1}(\bm{\theta}_{i-1})|\mathcal{D}_{1}^{N-1}\right]$
$\displaystyle\leq\frac{4}{N}\inf_{f\in\mathcal{F}(\bm{\theta}_{0},R)}\sum_{i=1}^{N}\psi(y_{i}f_{i})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}+\mathcal{O}\Big{(}\sqrt{\frac{\log(1/\delta)}{N}}\Big{)}+\Delta$

where $\widetilde{\bm{\theta}}$ is uniformly sampled from
$\\{\bm{\theta}_{1},\ldots,\bm{\theta}_{N-1}\\}$. ∎

###  C.3 Proof of Theorem 2

In the following, we show that the expected error is bounded by
$\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}}$ and is
proportional to $L((3\rho+1)\tau)^{L-1}/\sqrt{N}$.

######  Theorem 2 (Multi-layer GNN-based method).

For any $\delta\in(0,1/e]$ and $R>0$, there exists there exists
$m^{\star}=\mathcal{O}\left(N^{2}/L^{2}\right)\log(1/\delta),$ such that if
$m\geq m^{\star}$, then with probability at least $1-\delta$ over the
randomness of $\bm{\theta}_{0}$ with step size
$\eta=\frac{R}{m\sqrt{2N}((3\rho+1)\tau)^{(L-1)}}$ we have
$\mathbb{E}\left[\text{loss}_{N}^{0-1}(\widetilde{\bm{\theta}})|\mathcal{D}_{1}^{N-1}\right]\leq\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}+\mathcal{O}\Big{(}\sqrt{\frac{\log(1/\delta)}{N}}\Big{)}+\Delta,$
where
$R=\mathcal{O}(\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}})$,
$\widetilde{\bm{\theta}}$ is uniformly selected from
$\\{\bm{\theta}_{0},\ldots,\bm{\theta}_{N-1}\\}$, $\mathcal{D}_{1}^{N-1}$ is
the sequence of data points sampled before the $N$-th iteration, and the
expectation is computed on the uniform selection of weight parameters
$\widetilde{\bm{\theta}}$ and the condition sampling of $N$-th iteration data
examples.

The proof of Theorem 2 follows the proof of Corollary 3.10 in [Cao and Gu,
2019].

To begin with, let us recall from Lemma 6 that

$\displaystyle\mathbb{E}\left[\text{loss}_{N}^{0-1}(\bm{\theta})|\mathcal{D}_{1}^{N-1}\right]$
$\displaystyle\leq\frac{4}{N}\inf_{f\in\mathcal{F}(\bm{\theta}_{0},R)}\sum_{i=1}^{N}\psi(y_{i}f_{i})+\mathcal{O}\Big{(}\frac{LR((3\rho+1)\tau)^{(L-1)}}{\sqrt{N}}\Big{)}+\mathcal{O}\left(\sqrt{\frac{\log(1/\delta)}{N}}\right)+\Delta.$

Our goal is to show that $\exists f^{\prime}\in\mathcal{F}(\bm{\theta}_{0},R)$
such that
$\psi(y_{i}f^{\prime}_{i})=\log(1+\exp(-y_{i}f_{i}))\leq\frac{1}{\sqrt{N}},\leavevmode\nobreak\
\forall v_{i}\in\mathcal{V}$, which implies

$y_{i}f^{\prime}_{i}\geq-\log(\exp(N^{-1/2})-1),\leavevmode\nobreak\ \forall
v_{i}\in\mathcal{V}.$

Let us define $B=-\log(\exp(N^{-1/2})-1)$ and
$B^{\prime}=\max_{i\in\mathcal{V}}|f_{i}(\bm{\theta}_{0})|$ for notation
simplicity.

By the definition of $\mathcal{F}(\bm{\theta}_{0},R)$, we know that
$\exists\leavevmode\nobreak\
\bm{\theta}\in\mathcal{B}(\bm{\theta}_{0},Rm^{-1/2})$ such that

$\displaystyle
y_{i}f_{i}^{\prime}=y_{i}\Big{(}f_{i}(\bm{\theta}_{0})+\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}\rangle\Big{)}$
$\displaystyle=y_{i}f_{i}(\bm{\theta}_{0})+y_{i}\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}\rangle$
$\displaystyle\geq-B^{\prime}+y_{i}\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}\rangle\geq
B,$

where the inequality holds because

$-\max_{i\in\mathcal{V}}|f_{i}(\bm{\theta}_{0})|\leq
y_{i}f_{i}(\bm{\theta}_{0})\leq\max_{i\in\mathcal{V}}|f_{i}(\bm{\theta}_{0})|\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\
\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}\rangle=y_{i}(B+B^{\prime}).$

For notation simplicity, we define

$\displaystyle\hat{y}_{i}$
$\displaystyle=\langle\nabla_{\theta}f_{i}(\bm{\theta}_{0}),\bm{\theta}\rangle,\leavevmode\nobreak\
$ $\displaystyle\hat{y}_{i}$
$\displaystyle=(B+B^{\prime})y_{i},\leavevmode\nobreak\ $ $\displaystyle M$
$\displaystyle=md+m+m^{2}(L-2).$

Besides, let us denote the stack of gradient is

$\mathbf{J}=\begin{bmatrix}\text{vec}(\nabla_{\theta}f_{1}(\bm{\theta}_{0}))\\\
\text{vec}(\nabla_{\theta}f_{1}(\bm{\theta}_{1}))\\\ \vdots\\\
\text{vec}(\nabla_{\theta}f_{1}(\bm{\theta}_{N}))\end{bmatrix}\in\mathbb{R}^{N\times
M},\leavevmode\nobreak\ M=|\bm{\theta}|$

Besides, let us define
$\mathbf{J}=\mathbf{P}\mathbf{\Lambda}\mathbf{Q}^{\top}$ as the singular value
decomposition of $\mathbf{J}$, where $\mathbf{P}\in\mathbb{R}^{N\times
N},\mathbf{Q}\in\mathbb{R}^{M\times M}$ have orthonormal columns, and
$\mathbf{\Lambda}\in\mathbb{R}^{N\times M}$ is the singular value matrix.

Let us define
$\mathbf{w}=\mathbf{Q}\mathbf{\Lambda}^{-1}\mathbf{P}^{\top}\hat{\mathbf{y}}$,
then multiplying both sides by $\mathbf{J}$ we have

$\mathbf{J}\mathbf{w}=(\mathbf{P}\mathbf{\Lambda}\mathbf{Q}^{\top})(\mathbf{Q}\mathbf{\Lambda}^{-1}\mathbf{P}^{\top})\hat{\mathbf{y}}=\hat{\mathbf{y}}.$

Since
$\|\hat{\mathbf{y}}\|_{2}^{2}=\|(B+B^{\prime})\mathbf{y}\|_{2}^{2}=(B+B^{\prime})^{2}N$,
we have

$\displaystyle\|\mathbf{w}\|_{2}^{2}$
$\displaystyle=\|\mathbf{Q}\mathbf{\Lambda}^{-1}\mathbf{P}^{\top}\hat{\mathbf{y}}\|_{2}^{2}$
$\displaystyle=\hat{\mathbf{y}}^{\top}\mathbf{P}\mathbf{\Lambda}^{-1}\mathbf{Q}^{\top}\mathbf{Q}\mathbf{\Lambda}^{-1}\mathbf{P}^{\top}\hat{\mathbf{y}}$
$\displaystyle=\hat{\mathbf{y}}^{\top}\mathbf{P}\mathbf{\Lambda}^{-2}\mathbf{P}^{\top}\hat{\mathbf{y}}$
$\displaystyle=\hat{\mathbf{y}}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\hat{\mathbf{y}}$
$\displaystyle=(B+B^{\prime})^{2}N\cdot\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}.$

Let $\bm{\theta}$ be the parameters reshaped from $\mathbf{w}$, then we have

$\|\bm{\theta}\|_{\mathrm{F}}\leq\|\mathbf{w}\|_{2}\leq\mathcal{O}\Big{(}\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}}\Big{)}\Rightarrow\bm{\theta}\in\mathcal{B}\Big{(}\bm{0},\mathcal{O}\Big{(}\sqrt{\mathbf{y}^{\top}(\mathbf{J}\mathbf{J}^{\top})^{-1}\mathbf{y}}\Big{)}\Big{)},$

which concludes our proof.

##  Appendix D Generalization bound of RNN-based method

Recall that we compute the representation of node $v_{i}$ at time $t$ by
applying a multi-step RNN onto a sequence of temporal events
$\\{\mathbf{v}_{1},\ldots,\mathbf{v}_{L-1}\\}$ that are constructed at the
target node. The temporal event features $\mathbf{v}_{\ell}$ are pre-computed
on the temporal graph.

Representation computation. The RNN has trainable parameters
$\bm{\theta}=\\{\mathbf{W}^{(1)},\mathbf{W}^{(2)},\mathbf{W}^{(3)}\\}$ and
$\alpha\in\\{0,1\\}$ is binary hyper-parameter that controls whether a
residual connection is used. Then, the final prediction on node $v_{i}$ is
computed as $f_{i}(\bm{\theta})=\mathbf{W}^{(4)}\mathbf{h}_{L-1}$, where the
hidden representation $\mathbf{h}_{L-1}$ is recursively compute by

$\mathbf{h}_{\ell}=\sigma\Big{(}\kappa\big{(}\mathbf{W}^{(1)}\mathbf{h}_{\ell-1}+\mathbf{W}^{(2)}\mathbf{x}_{\ell}\big{)}\Big{)}+\alpha\mathbf{h}_{\ell-1}\in\mathbb{R}^{m}.$

Here $\sigma(\cdot)$ is the activation function,
$\mathbf{h}_{0}=\mathbf{0}_{m}$ is initialized as all-zero vector, and
$\mathbf{x}_{\ell}=\mathbf{W}^{(0)}\mathbf{v}_{\ell}$. We normalize the hidden
representation by $\kappa=1/\sqrt{2}$ so that $\|\mathbf{h}_{\ell}\|_{2}^{2}$
does not grow exponentially with respect to the number of steps $L$. For
weight parameters, we have
$\mathbf{W}^{(1)},\mathbf{W}^{(2)}\in\mathbb{R}^{m\times m}$,
$\mathbf{W}^{(3)}\in\mathbb{R}^{1\times m}$ as the trainable parameters, but
$\mathbf{W}^{(0)}\in\mathbb{R}^{m\times d}$ is non-trainable.

Gradient computation. The gradient with respect to each weight matrix is
computed by

$\displaystyle\frac{\partial f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(1)}}$
$\displaystyle=\sum_{\ell=1}^{L-1}\kappa\left[\mathbf{W}^{(3)}\left(\kappa\mathbf{D}_{L-1}\mathbf{W}^{(1)}+\alpha\mathbf{I}_{m}\right)\ldots\left(\kappa\mathbf{D}_{\ell+1}\mathbf{W}^{(1)}+\alpha\mathbf{I}_{m}\right)\mathbf{D}_{\ell}\right]^{\top}\mathbf{h}_{\ell}^{\top}\in\mathbb{R}^{m\times
m}$ $\displaystyle\frac{\partial
f_{i}(\bm{\theta})}{\partial\mathbf{W}^{(2)}}$