Efficient and Effective Implicit Dynamic Graph Neural Network

Yongjian Zhong Hieu Vu

Department of Computer Science Department of Computer Science
University of Iowa University of Iowa
Iowa City, IA, USA Iowa City, IA, USA
yongjian-zhong@[Link] hieu-vu@[Link]

Tianbao Yang Bijaya Adhikari

Department of Computer Science and Engineering Department of Computer Science
Texas A&M University University of Iowa
arXiv:2406.17894v1 [[Link]] 25 Jun 2024

College Station, TX, USA Iowa City, IA, USA

tianbao-yang@[Link] bijaya-adhikari@[Link]

ABSTRACT ACM Reference Format:

Implicit graph neural networks have gained popularity in recent Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari. 2024. Efficient
and Effective Implicit Dynamic Graph Neural Network. In Proceedings of
years as they capture long-range dependencies while improving
ACM Conference (Conference’17). ACM, New York, NY, USA, 12 pages. https:
predictive performance in static graphs. Despite the tussle between //[Link]/10.1145/[Link]
performance degradation due to the oversmoothing of learned em-
beddings and long-range dependency being more pronounced in
dynamic graphs, as features are aggregated both across neighbor-
1 INTRODUCTION
hood and time, no prior work has proposed an implicit graph neural Graph Convolution Network (GCN) [12] and its subsequent vari-
model in a dynamic setting. ants [16] [32] have achieved the state-of-the-art performance in pre-
In this paper, we present Implicit Dynamic Graph Neural Net- dictive tasks in various applications including molecular prediction
work (IDGNN) a novel implicit neural network for dynamic graphs [23], recommendation [17], and hyperspectral image classification
which is the first of its kind. A key characteristic of IDGNN is that [9]. GCNs have also been extended to the dynamic setting, where
it demonstrably is well-posed, i.e., it is theoretically guaranteed to the graph changes over time. Even in the dynamic setting, GCNs
have a fixed-point representation. We then demonstrate that the have achieved state-of-the-art results for numerous tasks including
standard iterative algorithm often used to train implicit models is rumor detection [31] and traffic prediction [13].
computationally expensive in our dynamic setting as it involves Despite numerous advantages, a major limitation of existing
computing gradients, which themselves have to be estimated in an GCNs is that they can only aggregate information up to 𝑘-hops,
iterative manner. To overcome this, we pose an equivalent bilevel where 𝑘 is the depth of the graph convolution operation. Hence,
optimization problem and propose an efficient single-loop training standard graph neural networks [12, 32] cannot capture long-range
algorithm that avoids iterative computation by maintaining moving dependencies beyond a radius imposed by the number of convo-
averages of key components of the gradients. We conduct exten- lution operations used. Trivial solutions like setting 𝑘 to a large
sive experiments on real-world datasets on both classification and number fail in overcoming this issue as empirical evidence [15]
regression tasks to demonstrate the superiority of our approach suggests that deepening the layers of GCN, even beyond a few
over the state-of-the-art baselines. We also demonstrate that our (2-4) layers, can lead to a notable decline in their performance.
bi-level optimization framework maintains the performance of the This is because the stacked GCN layers gradually smooth out the
expensive iterative algorithm while obtaining up to 1600x speed- node-level features which eventually results in non-discriminative
up. embeddings (aka oversmoothing). This creates a dilemma where,
on the one hand, we would like to capture dependencies between
nodes that are far away in the network by stacking multiple layers
KEYWORDS
of GCN together. On the other hand, we also would like to maintain
Dynamic Graphs, Implicit Graph Neural Networks, Graph Convo- the predictive performance by only using a few layers. To tackle
lutional Networks this dilemma in the static setting, Gu et al. [7] proposed an implicit
graph neural network (IGNN), which iterates the graph convolution
operator until the learned node representations converge to a fixed-
Permission to make digital or hard copies of all or part of this work for personal or
classroom use is granted without fee provided that copies are not made or distributed point representation. Since there is no a priori limitation on the
for profit or commercial advantage and that copies bear this notice and the full citation number of iterations, the fixed-point representation potentially con-
on the first page. Copyrights for components of this work owned by others than the tains information from all neighbors in the graph. Evidence shows
author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or
republish, to post on servers or to redistribute to lists, requires prior specific permission that it is able to capture long-range dependency while maintaining
and/or a fee. Request permissions from permissions@[Link]. predictive performance. Following this, other recent works [18, 22]
Conference’17, July 2017, Washington, DC, USA also have addressed the problem in the static setting.
© 2024 Copyright held by the owner/author(s). Publication rights licensed to ACM.
ACM ISBN 978-x-xxxx-xxxx-x/YY/MM In the case of dynamic graphs, where graphs evolve over time,
[Link] dynamic graph neural networks aggregate information over the
Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

current graph topology and historical graphs to learn meaningful components: propagation and update, which enable information
representations [31, 36]. Note the architecture of the graph neural aggregation and propagation for new interactions. EvolveGCN [21]
networks within each time stamp in dynamic graph neural net- uses an RNN to update GCN parameters and capture dynamic graph
works is similar to existing static GCNs. Hence, the information properties. Sankar et. al. [29] propose a Dynamic Self-Attention
aggregated for each node in a given time stamp is still limited to a Network (DySAT) with structural and temporal blocks to capture
radius of 𝑘 −ℎ𝑜𝑝𝑠 within the time stamp. Increasing the depth of the graph information. TGN [26] models edge streaming to learn node
GCN operator in each time stamp to ensure long-range dependency embeddings using an LSTM for event memory. TGAT [35] consid-
exacerbates the performance degradation of dynamic graph neural ers the time ordering of node neighbors. Gao et al. [5] explores the
networks as these models convolve both over the time stamps and expressiveness of temporal GNN models and introduces a time-then-
within the time stamps, hence oversmoothing the features even graph framework for dynamic graph learning, leveraging expressive
faster. Therefore, capturing long-range dependency while improv- sequence representations like RNN and transformers.
ing (or even maintaining) the performance is a big challenge for Implicit Graph Models: The implicit models or deep equilib-
dynamic graph neural networks. Despite its importance, very few rium models define their output using fixed-point equations. [1].
prior works have studied this phenomenon in dynamic graphs: Yang propose an equilibrium model for sequence data based on the fixed-
et al. [36] propose an L2 feature normalization process to alleviate point solution of an equilibrium equation. El et al. [4] introduce
the smoothing of features in dynamic graphs and Wang et al. [33] a general implicit deep learning framework and discuss the well-
mitigate the oversmoothing problem by emphasizing the impor- posedness of implicit models. Gu et al. [7] demonstrate the potential
tance of low-order neighbors via a node-wise encoder. However, of implicit models in graph representation learning, specifically
these approaches either rescale features or forget neighborhood with their implicit model called IGNN, which leverages a few layers
information, both of which are not ideal. of graph convolution network (GCN) to discover long-range de-
To address the challenges mentioned above, we propose IDGNN, pendencies. Park et al. [22] introduce the equilibrium GNN-based
an implicit neural network for dynamic graphs derived from the model with a linear transition map, and they ensure the transition
first principles. In designing IDGNN, we encountered multiple chal- map is contracting such that the fixed point exists and is unique.
lenges including i) uncertainty on whether fixed-point (converged) Liu et al. [18] propose an infinite-depth GNN that captures long-
representations exist for implicit neural models defined over dy- range dependencies in the graph while avoiding iterative solvers
namic graphs; and ii) efficiently training a model to find these by deriving a closed-form solution. Chen et al. [2] employ the dif-
fixed-point representations. In this paper, we overcome the first fusion equation as the equilibrium equation and solve a convex
challenge by providing theoretical guarantees on the existence optimization problem to find the fixed point in their model.
of the fixed-point representations on a single dynamic graph by Implicit Models Training: Efficiently training implicit models
leveraging a periodic model and generalizing this result to a set has always been a key challenge. Normally, the gradient of im-
of dynamic graphs. For the second challenge, we notice that the plicit models is obtained by solving an equilibrium equation using
stochastic gradient descent via implicit differentiation (often used fixed-point iteration or reversing the Jacobian matrix [7]. However,
by other implicit models [7, 14]) is too inefficient in our problem set- training these models via implicit deferential introduces additional
ting. As such, we reformulate our problem as an equivalent bilevel computational overhead. Geng et al. [6] propose a phantom gra-
optimization problem and design an efficient optimization strategy. dient to accelerate the training of implicit models based on the
The key contributions of the paper are as follows: damped unrolling and Neumann series. Li et al. [14] leverage sto-
chastic proximal gradient descent and its variance-reduced version
• We propose a novel dynamic graph neural network IDGNN,
to accelerate the training.
which ensures long-range dependency while providing theo-
retical guarantees on the existence of fixed-point representa-
tions. IDGNN is the first approach to leverage implicit graph
neural network framework for dynamic graphs. 3 METHODOLOGY
• We present a bilevel optimization formulation of our problem
and propose a novel stochastic optimization algorithm to 3.1 Preliminaries
efficiently train our model. Our experiments show that the Dynamic Graphs: We are given a set of 𝑁 dynamic graphs {G𝑖 }𝑖=1 𝑁 .
proposed optimization algorithm is faster than the naive 1 𝑡 𝑇
Each dynamic graph G𝑖 = {𝐺𝑖 , ..., 𝐺𝑖 , ..., 𝐺𝑖 } is a collection of 𝑇
gradient descent by up to 1600 times. snapshots. Let V denote the union set of nodes that appear in
• We conduct comprehensive comparisons with existing meth- any dynamic graph and 𝑛 := |V | be the total number of nodes.
ods to demonstrate that our method captures the long-range Without loss of generality, each snapshot can be represented as
dependency and outperforms the state-of-the-art dynamic 𝐺𝑖𝑡 = {V, E𝑖𝑡 , 𝑋𝑖𝑡 } since we can assume each graph is built on the
graph neural models on both classification and regression union set V and treat the absent nodes as isolated. E𝑖𝑡 is the set
tasks. of edges at time 𝑡 in G. 𝑋𝑖𝑡 ∈ R𝑙 ×𝑛 represents the node attribute
matrix, where 𝑙 is the dimension of the node attributes. Let 𝐴𝑖𝑡 be
2 RELATED WORK the adjacency matrix of 𝐺𝑖𝑡 .
Dynamic Graph Representation Learning: GNN has been suc- In this paper, we focus on node-level tasks (e.g. classification
cessful for static graphs, leading to the development of GNN-based and regression) for dynamic graphs. We consider the dataset as
algorithms for dynamic graphs [11]. DyGNN [19] comprises two 𝑁 , where G is the 𝑖-th dynamic graph, and 𝒚 ∈ R𝑛 is
{(G𝑖 , 𝒚𝑖 )}𝑖=1 𝑖 𝑖
Efficient and Effective Implicit Dynamic Graph Neural Network Conference’17, July 2017, Washington, DC, USA

for a shared 𝑉 for simplicity (thorough empirical discussion on this

choice is presented in the Experiment section). Following the prin-
ciple of the implicit model [1, 4, 7], we apply our model iteratively
infinite times. If the process converges, we consider the converged
result {𝑍 ∞1 , . . . , 𝑍 𝑇 } as the final embeddings. Consequently, the
∞
final embeddings have to satisfy the system of equations in (1) and
can be considered a fixed-point solution to (1). However, at this
point, it is not clear whether the fixed-point solution always exists
for arbitrary graph G.
Well-posedness is a property that an implicit function, such as
in (1), possesses a unique fixed point solution. While Gu et al. [7]
demonstrated the well-posedness of a single-layer implicit GCN
on a single static graph, the question remains open for dynamic
graphs. To establish the well-posedness property for our model, we
first introduce its vectorized version as follows.
1 𝜏 (1)
𝑧𝑘+1 = 𝜎 (𝑀 1𝑧𝑘 + vec(𝑉 𝑋 1 ))
2 𝜏 (2)
Figure 1: Model overview. This figure indicates the forward 𝑧𝑘+1 = 𝜎 (𝑀 2𝑧𝑘 + vec(𝑉 𝑋 2 ))
process of our model. ···
𝜏 (𝑇 )
𝑧𝑇𝑘+1 = 𝜎 (𝑀𝑇 𝑧𝑘 + vec(𝑉 𝑋 𝑇 )) (2)
the node-level labels assigned to the nodes in the last snapshot of
dynamic graph G𝑖 . where 𝑧 = vec(𝑍 ) is column-wise vectorization of 𝑍 , and = 𝑀𝑖
Permutation: Here, we define a permutation function to simplify (𝐴𝑖 ) ⊤ ⊗𝑊 𝑖 where ⊗ is the Kronecker product. Note that Equations
the notation in the rest of the paper. Let 𝜏 (𝑡) := 𝑇 − [(𝑇 − 𝑡 + 1) (2) can also be expressed in a single matrix form. This transforma-
mod 𝑇 ], where 𝑇 refers to the number of snapshots. This function tion involves sequentially connecting the shared nodes between
maps 𝑡 → 𝑡 − 1 for 𝑡 ∈ [2,𝑇 ] and maps 𝑡 → 𝑇 for 𝑡 = 1. the graphs. Thus, the formula (2) can be reformulated as follows:
Implicit Models: In general, implicit models [4, 6, 7] have the form
𝑍 = 𝑓 (𝑍, 𝑋 ), where 𝑓 is a function that can be parameterized by a 𝑧1   0 0 ··· 0 𝑀 1  𝑧 1   vec(𝑉 𝑋 1 ) 
 𝑧   vec(𝑉 𝑋 2 )  ª®
 2 © 𝑀 2   2  
neural network, 𝑋 is the data and 𝑍 is the learned representation. 𝑧  ­ 0 ··· 0 0 
 𝑧   vec(𝑉 𝑋 3 )  ®®
 3
𝑀3
 3 
We can obtain the fixed-point representation via iteration 𝑍 ∗ =

  = 𝜎 ­­  0
𝑧  ­ ··· 0 0   +
lim𝑘→∞ 𝑍𝑘+1 = lim𝑘→∞ 𝑓 (𝑍𝑘 , 𝑋 ) = 𝑓 (𝑍 ∗, 𝑋 ).
®
 .  ­  .. .. .. .. ..   .   .. ®
 ..  ­ . . . . .   ..   . ®
Thus, the key to designing an implicit model for dynamic graphs  
𝑧𝑇 
  
𝑧𝑇  vec(𝑉 𝑋 𝑇 ) 

is to provide a function 𝑓 that leads to converged representations.   « 0
 0 ··· 𝑀𝑇 0     ¬
(3)
3.2 Implicit Model for Dynamic Graphs Here, we omit the subscript for simplicity. Equation (3) represents
We first consider a single dynamic graph G = {𝐺 1, ..., 𝐺𝑇 } with 𝑇 a single equilibrium form of Equation (2). It can also be viewed as
snapshots and discover its well-posedness condition in this section. the time-expanded static version [] of our original dynamic graph
We then generalize the well-posedness conclusion for a set of dy- G. Based on the Banach fixed-point theorem [27], the Equation
namic graphs later. All the proofs are provided in the Appendix. (3) admits a unique fixed-point if the right-hand side is a contrac-
Now, let us consider the following stacked GCN model: tive mapping w.r.t. 𝑧. Therefore, we express the well-posedness
𝜏 (1) 1
condition for our model as follows,
1
𝑍𝑘+1 = 𝜎 (𝑊 1𝑍𝑘 𝐴 + 𝑉 𝑋 1)
𝜏 (2) 2
Theorem 3.1. For any element-wise non-expansive function 𝜎 (·),
2
𝑍𝑘+1 = 𝜎 (𝑊 2𝑍𝑘 𝐴 + 𝑉 𝑋 2) the coupled equilibrium equations in (2) have a unique fixed point
··· solution if ∥M ∥𝑜𝑝 < 1, where M define as
𝜏 (𝑇 ) 𝑇
𝑇
𝑍𝑘+1 = 𝜎 (𝑊 𝑇 𝑍𝑘 𝐴 + 𝑉 𝑋𝑇 ) (1)  0 ··· 0 𝑀 1
 2 
2 of the nodes
𝑀 ··· 0 0 
In the model presented above, the embeddings 𝑍𝑘+1 
 .. .. .. .. 
in the second time stamp in the (𝑘 + 1)-th layer depend on the  .
 . . . 
embeddings 𝑍𝑘1 of nodes in the first time stamp learned in the 𝑘-  0
 ··· 𝑀𝑇 0 
th layer and the feature of the nodes in the second time stamp and ∥M ∥𝑜𝑝 is the operator norm of M, which is the largest absolute
𝑋 2 . This design enables us to propagate information between time
eigenvalue. Furthermore, this is equivalent to ∥𝑀 𝑡 ∥𝑜𝑝 < 1 for any
stamps when stacking layers. The parameters for the 𝑡-th layer of
𝑡 = 1, ...,𝑇 .
the model are denoted as 𝑊 𝑡 ∈ R𝑑 ×𝑑 and 𝑉 ∈ R𝑑 ×𝑙 with 𝑉 being a
shared weight across all layers. Note that the proposed model and In order to maintain ∥𝑀 𝑡 ∥𝑜𝑝 < 1, ∀𝑡 = 1, ...,𝑇 , it is necessary to
the corresponding theory still hold when 𝑉 is not shared. We opt ensure that the condition 𝜆pr (|𝑊 𝑡 |)𝜆pr (𝐴𝑡 ) < 1 is satisfied, where
Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

𝜆pr (·) represents the Perron-Frobenius eigenvalue. However, satis- 4 TRAINING

fying the condition is challenging in general as it is hard to track In this section, we provide two algorithms to solve Equation (4). We
the eigenvalue as the underlying matrix changes. To overcome this first explore the stochastic gradient descent (SGD) method where
challenge, we impose a more stringent requirement on 𝑊 which we estimate the gradients leveraging the Implicit Function Theorem.
is more easily enforceable by leveraging a convex projection. We This approach offers several advantages, such as eliminating the
formally state this in the following theorem. need to store intermediate results during the forward pass and en-
Theorem 3.2. Let 𝜎 be an element-wise non-expansive function. abling direct backpropagation through the equilibrium point. While
If the coupled equilibrium equations satisfy the well-posedness con- widely used in various techniques [1, 7], this approach presents
dition, namely ∥M 𝑡 ∥𝑜𝑝 ≤ ∥𝑊 𝑡 ∥𝑜𝑝 ∥𝐴𝑡 ∥𝑜𝑝 < 1, ∀𝑡 = 1, ...,𝑇 , then certain drawbacks when applied to our specific model, particu-
there exists rescale coupled equilibrium equations, which satisfy the larly regarding computational overhead. Subsequently, we intro-
condition ∥𝑊 𝑡 ∥ ∞ ∥𝐴𝑡 ∥𝑜𝑝 < 1, ∀𝑡 = 1, ...,𝑇 , and the solutions of these duce an efficient training algorithm for our model, which adopts a
two equations are equivalent. bilevel viewpoint of our problem. This novel approach allows us to
overcome the limitations of the vanilla SGD method, resulting in
As previously stated, the theorem is formulated for a single improved computational efficiency during training.
dynamic graph. In the following Remark, we present a broader and
more general conclusion for a set of dynamic graphs. 4.1 SGD with Implicit Differentiation
Remark 1. Considering a set of dynamic graphs denoted as {G𝑖 }𝑖=1 𝑁 , The algorithm operates as follows: it first finds the fixed-point
achieving fixed-point representations for all dynamic graphs using embedding through iteration and then computes the gradient based
our model is guaranteed if the condition ∥𝑊 𝑡 ∥ ∞ ∥𝐴𝑖𝑡 ∥𝑜𝑝 < 1 holds on this embedding. It then updates the weights using SGD. The
for all time steps 𝑡 = 1, ...,𝑇 and all graph indices 𝑖 = 1, ..., 𝑁 . main obstacle in our way is estimating the gradient.
The parameters that need to be updated are 𝑊 , 𝑉 (from GCN
In practice, we can track the matrices with the largest operator
layers), and 𝜃 (from the classifier). The gradient with respect to
norms (max𝑖 ∥𝐴𝑖𝑡 ∥𝑜𝑝 ) at each time step. Focusing on these "critical"
parameter 𝜃 can be obtained as 𝜕𝜕𝜃L , which can be computed using
matrices and their corresponding constraints ensures our model
autograd functions given the fixed point. However, computing the
meets the required conditions. 𝜕L
gradient for other parameters presents a greater challenge. Let 𝜕𝑃
Incorporating Downstream Tasks: Based on the established con- 𝑖

ditions, we can obtain the fixed-point representation by iteratively represent the gradient with respect to 𝑊𝑖 or 𝑉𝑖 . Then, the gradient
𝜕L Í 𝜕 L 𝜕𝑧 𝑗 𝜕L
applying our model. Now, we want the fixed-point representations is computed as 𝜕𝑃 𝑖
= 𝑇𝑗=1 𝜕𝑧 𝑗 𝜕𝑃𝑖
. The computation of 𝜕𝑧 𝑗
can be
suited for specific downstream tasks with their own optimization achieved through the autograd mechanism. However, determining
𝜕𝑧 𝑗
objective. Let us now introduce the comprehensive objective which 𝜕𝑃𝑖 is non-trivial due to the cyclic definition of 𝑧 𝑗 .
incorporates both the application loss and the convergence require- By the definition of our model, the learned embeddings must
ments mentioned above. To this end, we utilize a neural network satisfy the following equations:
𝑓𝜃 (.), parameterized by 𝜃 , to map graph embeddings to their respec-
tive targets. Let W := {𝑊 1, ...,𝑊 𝑇 }. The comprehensive objective 𝑧 1 − 𝜎 (𝑀 1𝑧𝜏 (1) + vec(𝑉 𝑋 1 )) = 0
can now be summarized as follows: 𝑧 2 − 𝜎 (𝑀 2𝑧𝜏 (2) + vec(𝑉 𝑋 2 )) = 0
𝑁 ..
∑︁ .
min L (𝜃, W, 𝑉 ) = ℓ (𝑓𝜃 (𝑧𝑇𝑖 ), 𝒚𝑖 ) (4)
𝜃,W,𝑉
𝑖=1 𝑧𝑇 − 𝜎 (𝑀𝑇 𝑧𝜏 (𝑇 ) + vec(𝑉 𝑋 𝑇 )) = 0 (5)
  
𝜏 (1)
s.t. 𝑧𝑖1 =𝜎 (𝐴𝑖1 ) ⊤ ⊗𝑊 1
𝑧𝑖 + vec(𝑉 𝑋𝑖1 ) We apply column-wise vectorization on matrices 𝑊 and 𝑉 respec-
tively to obtain 𝑤 and 𝑣. We can then calculate the gradients 𝜕𝑤 𝜕𝑧
··· 𝜕𝑧
   and 𝜕𝑣 using the implicit function theorem. Here, we will show the
𝜏 (𝑇 )
𝑧𝑇𝑖 = 𝜎 (𝐴𝑇𝑖 ) ⊤ ⊗ 𝑊 𝑇 𝑧𝑖 + vec(𝑉 𝑋𝑖𝑇 ) , details of the derivation.
𝜅 Let 𝑍𝑖𝑡𝑗 , 𝐴𝑖𝑡 𝑗 , 𝑋𝑖𝑡𝑗 and 𝑊𝑖𝑡𝑗 denote the element at the 𝑖-th row and 𝑗-
∥𝑊 𝑡 ∥ ∞ ≤ , ∀𝑖 = 1, ..., 𝑁 , 𝑡 = 1, ...,𝑇
∥𝐴𝑖𝑡 ∥𝑜𝑝 th column of 𝑍 𝑡 , 𝐴𝑡 , 𝑋 𝑡 and 𝑊 𝑡 , respectively. Taking the derivative
Í Í 𝜏 (𝑡 )
of element-wise form of Equation (5), 𝑍𝑖𝑡𝑗 −𝜎 ( 𝑙 𝑛 𝑊𝑖𝑙𝑡 𝑍𝑙𝑛 𝐴𝑛𝑡 𝑗 +
Í 𝑡 𝑎
Where ℓ is a loss function ( e.g. cross entropy loss, mean square 𝑙 𝑉𝑖𝑙 𝑋𝑙 𝑗 ) = 0, with respect to 𝑊𝑏𝑐 we obtain
error), and 𝜅 is a positive number that is close to 1, which is added 𝜏 (𝑡 )
for numerical stability. 𝜕𝑍𝑖𝑡𝑗 ∑︁
𝜏 (𝑡 ) 𝑡
∑︁ ∑︁ 𝜕𝑍𝑙𝑛
𝑎 − Σ 𝑡 ©
𝑖𝑗 ­ 𝛿 𝛿 𝑍
𝑎𝑡 𝑏𝑖 𝑐𝑛 𝐴 𝑛𝑗 + 𝑊 𝑡 𝑡 ª
𝑖𝑙 𝜕𝑊 𝑎 𝐴𝑛 𝑗 ® = 0
To find the fixed-point representation (i.e. forward pass), we can 𝜕𝑊𝑏𝑐
use fixed-point iteration or other root-finding algorithms. Backprop- «𝑛 𝑙 𝑛 𝑏𝑐 ¬
agation requires storing intermediate results, which is infeasible where 𝛿𝑎𝑡 is the indicator function which equals 1 only when 𝑎 = 𝑡,
since there might be hundreds of iterations in discovering the rep- and ⊙ denotes element-wise multiplication. Let 𝜎 ′ (.) represent the
resentation. Thus, the key challenge in solving Equation (4) lies in derivative 𝜎 (.), and we define
𝜏 (𝑡 )
∑︁ ∑︁ ∑︁
determining how to perform backpropagation effectively, especially Σ𝑡𝑖 𝑗 := 𝜎 ′ ( 𝑊𝑖𝑙𝑡 𝑍𝑙𝑛 𝐴𝑛𝑡 𝑗 + 𝑉𝑖𝑙 𝑋𝑙𝑡𝑗 )
in the context of dynamic graphs. 𝑙 𝑛 𝑙
Efficient and Effective Implicit Dynamic Graph Neural Network Conference’17, July 2017, Washington, DC, USA

which means Σ𝑡 is a matrix and has a shape like 𝑍 𝑡 . Let 𝐻 𝑡 := Algorithm 1 Bilevel Optimization for IGDNN
(𝑍 𝜏 (𝑡 ) 𝐴𝑡 )𝑇 , and 𝐻 ·𝑐
𝑡 denotes the 𝑐-th column of 𝐻 𝑡 . Let 𝑒 be a col-
𝑏 Require: D = {(G𝑖 , 𝒚𝑖 )}𝑖=1 𝑁 ,𝜂 ,𝜂 ,𝛾
1 2
umn vector with value 1 at 𝑏-th entry and 0 elsewhere. Enumerating Ensure: 𝜔, 𝜃
𝑍𝑖𝑡𝑗 for all 𝑖, 𝑗, we have 1: Randomly initialize 𝑧 0𝑗 and 𝑣 0𝑗 for 𝑗 = 1, ..., 𝑁
2: for 𝑘 = 0, 1, ..., 𝐾 do
!
𝜕𝑧𝑡 𝜏 (𝑡 )
𝑡 𝜕𝑧
𝑎 − vec(Σ ) ⊙ 𝛿𝑎𝑡 𝑒𝑏 ⊗ 𝐻 ·𝑐 + 𝑀 𝜕𝑊 𝑎
𝑡 𝑡
=0 3: Sample(a batch data B
𝜕𝑊𝑏𝑐 𝑏𝑐 (𝐼 − 𝜂 1 )𝑧ˆ𝑘𝑗 + 𝜂 1𝜙 (𝑧ˆ𝑘𝑗 , 𝜔 𝑡 ; G𝑖 ) 𝑗 ∈ B
4: 𝑧ˆ𝑘+1 =
Therefore, for any time stamp 𝑡, the gradient of 𝑧𝑡 with respect to 𝑗
𝑧ˆ𝑡 o.w.
the 𝑎-th layer of GCN, 𝑤 𝑎 , can be expressed as: ( 𝑗
(𝐼 − 𝜂 ∇ 2 𝑔(𝑧ˆ𝑘 , 𝜔 𝑘 )) 𝑣ˆ𝑘 + 𝜂 ∇ ℓ (𝑧ˆ𝑘 , 𝜔 𝑘 ) 𝑗 ∈ B
2 𝑧𝑧 𝑗 𝑗 2 𝑧 𝑗 𝑗
𝜏 (𝑡 )
! 5: 𝑣ˆ𝑘+1 = 𝑡
𝜕𝑧𝑡 𝑡 𝜕𝑧
𝑗
𝑣ˆ 𝑗 o.w.
− Ξ ⊙ 𝛿𝑎𝑡 𝐻 ⊗ 𝐼 + 𝑀
𝑡 𝑡
=0 (6)
𝜕𝑤 𝑎 𝜕𝑤 𝑎 6: Update gradient estimator
1 ∑︁ h i
Where Ξ𝑡 is a matrix that has identical column vectors, and each col- Δ𝑘+1 = ∇𝜔 ℓ 𝑗 (𝑧ˆ𝑘𝑗 , 𝜔 𝑘 ) − ∇𝜔𝑧2
𝑔 𝑗 (𝑧ˆ𝑘𝑗 , 𝜔 𝑘 )𝑣ˆ𝑘𝑗
umn vector is the vec(Σ𝑡 ). Similarly, we can compute the gradient |B|
𝑗∈B
of 𝑍𝑖𝑡𝑗 w.r.t. 𝑉𝑏𝑐 .
7: 𝑚𝑘+1 = (1 −𝛾)𝑚𝑘 + 𝛾 Δ𝑘+1
𝜏 (𝑡 )
𝜕𝑍𝑖𝑡𝑗 𝜕𝑍 8: 𝜔 𝑘+1 = Π Ω 𝜔 𝑘 − 𝜂 0𝑚𝑘+1
− Σ𝑡𝑖 𝑗 ­𝛿𝑏𝑖 𝑋𝑐𝑡 𝑗 + 𝑊𝑖𝑙𝑡 𝑙𝑛𝑎 𝐴𝑛𝑡 𝑗 ® =0
© ª
𝜕𝑉𝑏𝑐 𝜕𝑊𝑏𝑐 9: end for
« ¬
Therefore,
Where 𝜙 (𝑧, W, 𝑉 ; G𝑖 ) = 𝜎 (𝑀𝑖𝑇 ...𝜎 (𝑀𝑖1𝑧+vec(𝑉 𝑋𝑖1 ))...+vec(𝑉 𝑋𝑖𝑇 )).
!
𝜕𝑧𝑡 𝜕𝑧𝜏 (𝑡 )
− Ξ𝑡 ⊙ (𝑋 𝑡 ) ⊤ ⊗ 𝐼 + 𝑀 𝑡 =0 (7) The main differences between these problems lie in the constraints.
𝜕𝑣 𝜕𝑣
Equation (8) introduces explicit constraints solely on the last snap-
While these equations provide a path for gradient computation, shot, leading to a multi-block bilevel optimization problem. This
it is essential to note that all gradients are interconnected within type of problem has been investigated recently by [25] and [10].
a system of equations. The resolution of such a system entails a [25] focus on top-K NDCG optimization, formulating it as a com-
substantial computational overhead. positional bilevel optimization with a multi-block structure. Their
Per-iteration Complexity of naive gradient descent: Equa- approach simplifies updates by sampling a single block batch in
tions (6) and (7) reveal that a set of equilibrium equations determines each iteration and only updating the sampled blocks. [10] employs
the gradients. Consequently, to compute the gradients, we need to a similar technique but addresses a broader range of multi-block
solve these equations using fixed-point iteration. Each layer neces- min-max bilevel problems.
sitates one round of fixed-point iteration, and in total, including However, these state-of-the-art bilevel optimization algorithms
𝑉 , we need to perform fixed-point iteration 𝑇 + 1 times. The ma- are designed to address problems with strongly convex lower-level
jor computational overhead arises from the multiplication of 𝑀 problems, which does not hold true for our problem. It is evident
with the derivatives, resulting in a complexity of 𝑂 ((𝑛𝑑) 2𝑑 2 ). Each that our lower-level problems in 8 are generally nonconvex with
fixed-point iteration involves 𝑇 instances of such computations. respect to 𝑧 since they involve highly nonlinear neural networks.
Consequently, the overall runtime for each update is 𝑂 (𝑇 2𝑛 2𝑑 4 ). Additionally, these methods utilize stochastic gradient descent on
Although the adjacency matrix is sparse, it only reduces the com- the lower level in each iteration, which may lead to potential extra
plexity to 𝑂 (𝑇 2𝑛𝑑 4 ). This limits the application of our model in computation in estimating the gradient. Nevertheless, it is crucial
large-scale dynamic graphs and hampers our ability to utilize large to note that the optimal solution to our lower-level problem cor-
embeddings. responds to the fixed point of Eq (2). Leveraging this insight, we
employ a fixed-point iteration to update the lower-level solution.
4.2 Efficient Update via Bilevel Optimization We propose a single loop algorithm (see Algorithm 1) with fixed-
To address the previously mentioned challenges, we turn to Bilevel point updates.
Optimization [3] as a potential solution. We reformulate the prob- To better illustrate our algorithm, let 𝜔 = {W, 𝑉 }, and 𝑔𝑖 (𝑧, 𝜔)
lem presented in Equation (4) as the following standard bilevel represents the 𝑖-th-block lower-level problem, defined as ∥𝑧 −
optimization problem. 𝜙 (𝑧, 𝜔; G𝑖 )∥ 22 , and let ℓ𝑖 (𝑧, 𝜔) := ℓ (𝑓𝜃 (𝑧), 𝑦𝑖 ). The key to solving
the bilevel optimization problem in (8) is to estimate the hyper-
𝑁
∑︁ gradient ∇L (𝜃, 𝜔) and backpropagate. The hypergradient of our
min L (𝜃, W, 𝑉 ) = ℓ (𝑓𝜃 (𝑧𝑇𝑖 , 𝒚𝑖 ) (8) objective in (8) with respect to 𝜔 as follows:
𝜃,W,𝑉
𝑖=1 𝑁
1 ∑︁
s.t. 𝑧𝑇𝑖 = arg min ∥𝑧 − 𝜙 (𝑧, W, 𝑉 ; G𝑖 )∥ 22 ∇L (𝜃, 𝜔) = ∇ℓ𝑖 (𝑧𝑇𝑖 , 𝜔)
𝑧 𝑁 𝑖=1
𝑡 𝜅
∥𝑊 ∥ ∞ ≤ , ∀𝑖 = 1, ..., 𝑁 , 𝑡 = 1, ...,𝑇 2
h
2
i −1
∥𝐴𝑖𝑡 ∥𝑜𝑝 − ∇𝜔𝑧 𝑔𝑖 (𝑧𝑇𝑖 , 𝜔) ∇𝑧𝑧 𝑔𝑖 (𝑧𝑇𝑖 , 𝜔) ∇𝑧 ℓ𝑖 (𝑧𝑇𝑖 , 𝜔)
Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

Table 1: Performance for classification task (ROCAUC) and regression task (MAPE (%)). Performances on Brain10, England-
COVID, PeMS04, and PeMS08 for baseline methods are taken from [5]. The best performance for each dataset is highlighted in
bold, while the second-best performance is underlined.

Classification Regression
Model Brain10 DBLP5 Reddit4 England-COVID PeMS04 PeMS08
Trans. Induc. Trans. Induc. Trans. Induc.
EvolveGCN-O 0.58±0.10 0.639±0.207 0.513±0.008 4.07±0.73 3.88±0.47 3.20±0.25 2.61±0.42 2.65±0.12 2.40±0.27
EvolveGCN-H 0.60±0.11 0.510±0.013 0.508±0.008 4.14±1.14 3.50±0.42 3.34±0.14 2.84±0.31 2.81±0.28 2.81±0.23
GCN-GRU 0.87±0.07 0.878±0.017 0.513±0.010 3.56±0.26 2.97±0.34 1.60±0.14 1.28±0.04 1.40±0.26 1.07±0.03
DySAT-H 0.77±0.07 0.917±0.007 0.508±0.003 3.67±0.15 3.32±0.76 1.86±0.08 1.58±0.08 1.49±0.08 1.34±0.03
GCRN-M2 0.77±0.04 0.894±0.009 0.546±0.020 3.85±0.39 3.37±0.27 1.70±0.20 1.20±0.06 1.30±0.17 1.07±0.03
DCRNN 0.84±0.02 0.904±0.013 0.535±0.007 3.58±0.53 3.09±0.24 1.67±0.19 1.27±0.06 1.32±0.19 1.07±0.03
TGAT 0.80±0.03 0.895±0.003 0.510±0.011 5.44±0.46 5.13±0.26 3.11±0.50 2.25±0.27 2.66±0.27 2.34±0.19
TGN 0.91±0.03 0.887±0.004 0.521±0.010 4.15±0.81 3.17±0.23 1.79±0.21 1.19±0.07 1.49±0.26 0.99±0.06
GRU-GCN 0.91±0.03 0.906±0.008 0.525±0.006 3.41±0.28 2.87±0.19 1.61±0.35 1.13±0.05 1.27±0.21 0.89±0.07
IDGNN 0.94±0.01 0.907±0.005 0.556±0.017 2.65±0.25 3.05±0.25 0.53±0.05 0.63±0.04 0.45±0.11 0.50±0.05

 2  −1
Explicitly computing the inverted Hessian matrix ∇𝑧𝑧 𝑔𝑖 (𝑧𝑇𝑖 , 𝜔) IDGNN against nine state-of-the-art baselines on multiple real-
in computation of our hypergradient is computationally expen- world datasets, comprising three node classification and four node
sive. Inspired by [10] and [25], we instead directly approximate regression datasets. Key dataset statistics are summarized in Table
 2  −1
∇𝑧𝑧 𝑔𝑖 (𝑧𝑇𝑖 , 𝜔) ∇𝑧 ℓ𝑖 (𝑧𝑇𝑖 , 𝜔) using 𝑣𝑖 for each block by moving 2, with further details available in section 5.1. Due to the space
average estimation (line 5). More specifically, we track the optimal constraints, specifics on experimental setup and hyperparameters
point of min𝑣 12 𝑣 ⊤ ∇𝑧𝑧 2 𝑔 (𝑧𝑇 , 𝜔)𝑣 − 𝑣 ⊤ ∇ ℓ (𝑧𝑇 , 𝜔) for each block by
𝑖 𝑖 𝑧 𝑖 𝑖 are provided in Appendices A.1 and A.2, respectively. Our code is
maintaining 𝑣𝑖 . Let 𝑧ˆ𝑖 be a moving average approximation to the available for reproducibility. 1 .
optimal lower-level solution 𝑧𝑇𝑖 . We use a single fixed-point update
for 𝑧ˆ𝑖 (line 4). We do not want to update all blocks in every iteration Table 2: Statistics of datasets. 𝑁 : number of dynamic graphs,
since this is impractical when the number of blocks is large. To |𝑉 |: number of nodes, min |𝐸|: minimum number of edges,
address this issue, we use stochastic training. For sampled blocks, max |𝐸|: maximum number of edges,𝑇 : window size, 𝑑: feature
we update their 𝑧ˆ and 𝑣, ˆ and we compute the hypergradient (line 6). dimension, 𝑦 label dimension
Note that the multiplication ∇𝜔𝑧 2 𝑔 (𝑧ˆ𝑘 , 𝜔 𝑘 ) 𝑣ˆ𝑘 (in line 6) can also
𝑗 𝑗 𝑗
be efficiently computed using Hessian vector product [24]. .As a 𝑁 |𝑉 | min |𝐸 | max |𝐸 | 𝑇 𝑑 𝑦
result, the training time for our algorithm is proportional to normal Brain10 1 5000 154094 167944 12 20 10
backpropagation, eliminating the need for fixed-point iterations. DBLP5 1 6606 2912 5002 10 100 5
Note that in cases where the lower-level problem is strongly Reddit4 1 8291 12886 56098 10 20 4
convex, the errors introduced by these approximations are well- PeMS04 16980 307 680 680 12 5 3
contained [10]. Our lower-level problem admits a unique fixed PeMS08 17844 170 548 548 12 5 3
point, hence, employing fixed-point iteration becomes an efficient English-COVID 54 129 836 2158 7 1 1
means of attaining the optimal lower-level solution, akin to the
effectiveness of gradient descent under strong convexity. Therefore,
it is justifiable to assert that our approximations are effective in this 5.1 Datasets
scenario, with empirical evidence robustly endorsing their practical
Brain dataset is derived from a real-world fMRI brain scans dataset
efficacy. 2 . In this dataset, nodes represent brain tissues, and edges capture
Per-iteration Complexity of bilevel optimization: The main
nodes’ activation time similarity in the examined period. The node
computational overheads are updating 𝑣 and estimating gradient.
attributes are generated from the fMRI signals [34]. Given the tem-
Both steps are involved with estimating a huge Hassian matrix,
poral graph, our goal is to predict the functionality of each node,
but, in practice, we use Hessian vector product to avoid explicitly
which is one of the 10 categories.
computing the Hessian matrix. Therefore, the dominant runtime
DBLP is an extracted co-author network from DBLP website 3 ,
of bilevel optimization is three times backpropagation. Each back-
where nodes are authors and edges represent co-authorship rela-
propagation takes 𝑂 (𝑇𝑛𝑑 2 + 𝑇𝑛 2𝑑).
tionships. The extracted authors are from 5 research areas [34],
serving as our class labels. Note attributes are word2vec represen-
5 EXPERIMENTS tations of titles and abstracts from related papers.
In this section, we present the performance of IDGNN in various 1 [Link]
tasks, including effectiveness in capturing long-range dependencies 2 [Link]

and avoiding oversmoothing on a synthetic dataset. We benchmark 3 [Link]

Efficient and Effective Implicit Dynamic Graph Neural Network Conference’17, July 2017, Washington, DC, USA

Figure 2: The left and middle are accuracy and loss curves when using 10 layers. The x-axis is epochs, and the y-axis is accuracy
and cross entropy loss, respectively. The right plot represents the accuracy results of all baselines under different layer settings.

Reddit dataset is extracted from a popular online platform Red- settings, with the exception of the inductive case in England-COVID.
dit 4 [8], where nodes represent posts, and edges are constructed Our method demonstrates a significant improvement for PeMS04
between nodes having similar keywords. The node attributes are and PeMS08, particularly in the transductive learning scenario. In
word2vec representations of the comments of a post [34], and node comparison to the second-best method, our proposed model reduces
labels correspond to one of the 4 communities or “subreddits" to the error by over 1% of MAPE, but our model on the inductive
which the post belongs. learning scenario does not enjoy such improvement. We attribute
PeMS04 & PeMS08 These two datasets represent traffic sensor net- this difference to our model’s tendency to separate nodes, even
works in two distinct districts of California during various months. when they have the same labels and topology. We delve into this
In this context, each node symbolizes a road sensor, while the edge phenomenon in the Appendix B.
weights denote the geographical distance between two sensors
(with uniform weights across all snapshots). Every sensor captures 5.3 Classification
average traffic statistics such as flow, occupancy, and speed over a We conducted classification experiments on Brain10, DBLP5, and
5-minute interval. The datasets utilized in this study are identical Reddit4 datasets. Since these datasets consist of only one dynamic
to those in [5], where we exclusively forecast attribute values for a graph, we focused on testing the transductive case. The evaluation
single future snapshot. was done using the Area under the ROC Curve (AUC) metric. The
England-COVID This temporal graph dataset is taken from [5], average prediction AUC values and their corresponding standard
which is created by aggregating information from the mobility log deviations are presented in Table 1. Our proposed model achieved
released by Facebook under Data For Good program 5 for England the top rank in 2 out of 3 datasets and was the second best in the
[20]. The nodes are cities, and the edges are transportation statistics remaining dataset. These results demonstrate that our model suc-
between them. Given a temporal graph of length 7, which represents cessfully captures the long-range dependencies within the dynamic
the past week, we need to predict the infection rate of the next day. graphs, as reflected in the learned embeddings.
5.2 Regression
5.4 Long-range Dependency and
For node-level tasks, there are two evaluation paradigms: transduc-
Oversmoothing
tive and inductive. Transductive evaluation allows the model to
access the nodes and edges in testing data during training (i.e., only Long-range Dependency: The toy example aims to test the ability
the labels are hidden), while inductive evaluation involves testing of all approaches to capture long-range dependencies. The toy data
the model on new nodes and edges. In simpler terms, transductive we constructed consists of {5, 10, 15} snapshots, with each snapshot
evaluation separates training and testing sets based on nodes, while being a clique of 10 nodes. Each node has 10 associated attributes.
inductive evaluation separates them based on time stamps (i.e., The task is to classify nodes at the last snapshot, where each node
models are trained on past time stamps and tested on the future). represents its own class (i.e., there are a total of 10 classes). The node
Here, we conduct experiments on both settings. attributes consist of randomly generated numbers, except for the
The datasets we used for the regression task are England-COVID, first snapshot, which uses the one-hot representation of the class.
PeMS04, and PeMS08. We use the mean average percentage error Successful classification of this dataset requires effective informa-
(MAPE) as our evaluation metric. The results are reported in Table tion aggregation starting in the initial time stamp, propagating the
1 with mean MAPE and standard deviation. Our proposed method class label information over time, and avoiding oversmoothing as
outperforms other methods in both transductive and inductive the information is propagated. In this dataset, there are no testing
nodes; all nodes are used for training.
4 [Link] The training results are presented on Figure 2. Our model is
5 [Link] compared with TGCN [37]. We also propose two more modified
Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

Table 3: Memory and runtime comparison results for all Table 4: Evaluating the smoothness of embeddings and the
methods on reddit4 and DBLP5 datasets. We report the mem- AUC on Brain10.
ory usage using MB and runtime using seconds per window
. Dirichlet Energy ↑ AUC ↑
Reddit4 DBLP5 Static (3 layers) 10.660 89.79
Mem. Time Mem. Time Static (4 layers) 8.781 88.67
EvolveGCN-O 42 0.0649±0.017 52 0.0672±0.014 Static (5 layers) 3.399 81.49
EvolveGCN-H 52 0.0904±0.020 82 0.0997±0.037 W/o looping 8.957 85.61
GCN-GRU 221 0.0733±0.012 200 0.1142±0.045 IDGNN 35.299 94.87
DySAT-H 181 0.1613±0.056 165 0.1343±0.012
GCRN-M2 322 0.4345±0.080 319 0.4934±0.076
Table 5: Runtime and performance comparison between SGD
DCRNN 223 0.1697±0.019 278 0.2121±0.039
and Bilevel Methods.
TGAT 793 0.0750±0.014 338 0.0770±0.015
TGN 450 0.0417±0.004 233 0.0454±0.012
GRU-GCN 4116 0.0199±0.008 580 0.0161±0.007 Runtime (s/win) Performance
IDGNN 89 0.0291±0.007 75 0.0302±0.002 SGD Bilevel SGD Bilevel
Brain10 624 0.390 94.7 94.9
PeMS04 0.72 0.049 0.63 0.63
baselines: IGNN-GRU and TIGNN, which are obtained by replacing PeMS08 0.29 0.046 0.56 0.50
the GCNs within GCN-GRU [30] and TGCN by IGNN. We ensure England-COVID 0.092 0.030 2.97 3.05
the comparison is fair by ensuring a similar number of parameters
are used, and we test all models on {5, 10, 15} layers. All meth-
ods are trained for a maximum of 2000 epochs, followed by the • W/o loop: This model shares the structure of IDGNN but
hyper-parameter selection approach described in the Appendix. As does not enforce the learned embeddings to be a fixed point.
shown in the figure, we explored the impact of varying the number We evaluate these methods on the Brain10 dataset, and the results
of layers on both baseline models and our proposed model. We are presented in Table 4. The "static" method exhibits oversmooth-
documented the resulting accuracies accordingly. Notably, TGCN ing as the layer count increases, and a similar issue is observed for
exhibits a discernible pattern of over-smoothing, evidenced by a per- the "W/o loop" method, despite stacking one layer for each snap-
formance decline with increasing layers. Conversely, IGNN-GRU shot. These findings strongly indicate that our model effectively
and TIGNN do not demonstrate such susceptibility. Additionally, addresses oversmoothing in dynamic graph learning by leveraging
we observed challenges in data fitting for both methods, whereas the implicit model structure.
our model consistently achieved 100% accuracy across different
layer configurations. This experiment underscores the robustness 5.5 Efficiency
of our architecture in fully unleashing the potential of implicit
We compare runtime and performance between SGD and bilevel
models in dynamic graph learning.
optimization algorithms. To this end, we compare Brain10, England-
Moreover, we perform an alternative experiment, wherein we
COVID, PeMS04, and PeMS08. The results are summarized on the
shift label information from snapshot 1 to snapshot 5. This adjust-
Table 5. The results are computed by averaging the runtime of a
ment ensures consistent difficulty in leveraging label information
whole epoch with the number of dynamic graphs 𝑁 .
across all models. Due to space constraints, we provide the detailed
These methods have similar performance, but the runtime results
results in the Appendix; however, the overall conclusion remains
show that the bilevel optimization algorithm is much faster than the
unaffected.
SGD. Especially, in the Brain10 dataset, bilevel algorithm achieves
Oversmoothing: In this part, we employ Dirichlet energy to as-
1600 times of speedup compared with SGD. Furthermore, we notice
sess the smoothness of the learned embeddings [28]. The Dirichlet
that the ratio of runtimes in PeMS04 and PeMS08 is 0.72
0.29 = 2.48, and
energy measures the mean distance between nodes and is defined
as follows: the squared ratio of their number of nodes is ( 307 2
170 ) = 3.26. This
√︄ confirms our complexity result for SGD, which is quadratic regard-
1 ∑︁ ∑︁
𝐷𝐸 (𝑍 ) = ∥𝑍𝑖 − 𝑍 𝑗 ∥ 2 ing the number of nodes. On the other hand, the bilevel method
|V | exhibits only linear dependency. We also present the memory usage
𝑖 ∈ V 𝑗 ∈ N𝑖
and runtimes of all methods on Reddit4 and DBLP5. The memory ef-
Here, N𝑖 denotes the neighbors of node 𝑖. Low Dirichlet energy ficiency of implicit models comes from the fact that implicit models
means smooth embedding. We compare our model with two differ- can use few parameters and do not need to store the intermediate
ent baselines: results. However, we need to store intermediate results and back-
• Static: We create a static graph by averaging the adjacency propagate for our bi-level method. Due to the simple RNN-free
matrices of all snapshots and using the feature matrix of architecture of our method, our approach is competitive in runtime
the last snapshot as its feature. Subsequently, we learn the and memory. We provide a memory and runtime comparison on
embeddings using a vanilla multi-layer GCN. DBLP5 and Reddit4. The results are summarized in Table 3.
Efficient and Effective Implicit Dynamic Graph Neural Network Conference’17, July 2017, Washington, DC, USA

𝑂 (𝑇𝑑 2 + 𝑇𝑛𝑑 2 + 𝑡 𝐸𝑡 𝑑)
Í
EvolveGCN-O layer for static information, while the share-𝑉 method employs
EvolveGCN-H 2 2 Í
𝑂 (𝑇𝑑 + 𝑇𝑛𝑑 + 𝑡 𝐸𝑡 𝑑) one GCN and 𝑇 linear layers. The not-share method achieves the
GCN-GRU 2 2 Í
𝑂 (𝑇𝑑 + 𝑇𝑛𝑑 + 𝑡 𝐸𝑡 𝑑) best result, although the improvement is negligible. However, it
𝑂 (𝑇𝑑 2 + 𝑇𝑛𝑑 + 𝑡 𝐸𝑡 𝑑 2 )
Í
DySAT increases the parameter size, resulting in significant computational
𝑂 (𝑇𝑑 2 + 𝑇𝑛𝑑 2 + 𝑡 𝐸𝑡 𝑑)
Í
GCRN-M2 overhead. Hence, we opt for the current configuration, as it delivers
DCRNN 2 2 Í
𝑂 (𝑇𝑑 + 𝑇𝑛𝑑 + 𝑡 𝐸𝑡 𝑑) good performance while minimizing parameter size, especially
𝑂 (𝑇𝑛𝑑 + 𝑡 𝐸𝑡 𝑑 2 )
Í
TGAT since the number of attributes may exceed the hidden dimension.
𝑂 (𝐸𝑎𝑔𝑔𝑑 + 𝑇𝑛𝑑 2 + 𝑡 𝐸𝑡 𝑑 2 )
Í
TGN We additionally evaluate an alternative baseline, denoted as "w/o
GRU-GCN 𝑂 (𝐸𝑎𝑔𝑔𝑑 + 𝑇𝑛𝑑 2 + 𝑇 𝐸𝑎𝑔𝑔𝑑 2 ) loop," wherein we eliminate the looping structure from IDGNN.
𝑂 (𝑇𝑛𝑑 2 + 𝑡 𝐸𝑡 𝑑)
Í
IDGNN The obtained results reveal that this model exhibits the lowest
Table 6: Summary of complexities of all models performance, underscoring the efficacy of our proposed approach.

6 CONCLUSIONS
In this paper, we propose a novel implicit graph neural network
5.6 Complexity for dynamic graphs. As far as we know, this is the first implicit
To contrast with the emprical runtime, here we present the theoret- model on dynamic graphs. We demonstrate that the implicit model
ical time complexities for our approach and the baselines. For an we proposed has the well-posedness characteristic. We proposed a
arbitrary temporal graph, we denote the total number of snapshots standard optimization algorithm using the Implicit Function Theo-
as 𝑇 , the total number of nodes as 𝑛, the number of edges of time 𝑡 rem. However, the optimization was too computationally expensive
as 𝐸𝑡 , and 𝐸𝑎𝑔𝑔 as the number of aggregated edges, which satisfies for our model. Hence, we proposed a novel bilevel optimization
Í Í
𝐸𝑎𝑔𝑔 ≪ 𝑡 𝐸𝑡 or 𝐸𝑎𝑔𝑔 ≈ 𝑡 𝐸𝑡 . Some basic complexities: GRU, self- algorithm to train our proposed model. We conducted extensive
attention, and LSTM have complexity 𝑂 (𝑇 2 ) if the input sequence experiments on 6 real-world datasets and one toy dataset. The re-
has length 𝑇 ; GCN and SpectralGCN have complexity 𝑂 (𝑛𝑑 2 + 𝐸𝑑); gression and classification tasks show that the proposed approach
GAT has complexity 𝑂 (𝑛𝑑 + 𝐸𝑑 2 ). The complexities of all models outperforms all the baselines in most settings. Finally, we also
are summarized in Table 6. Based on the complexity results, IDGNN demonstrated that the proposed bilevel optimization algorithm
is faster than EvolveGCN-O, EvolveGCN-H, GCN-GRU, GCRN-M2, obtains significant speedup over standard optimization while main-
and DCRNN due to the absence of RNN in our model. taining the same performance. A key limitation of our proposed
approach is that it can only predict the consecutive snapshot. In the
5.7 Ablation Study future, we plan on addressing this issue and also provide a diffusion
In this section, we will delve into the various configurations of model-based training algorithm.
our model. Drawing from the properties mentioned earlier, our
model can be represented by the formula (1). We have deliberately
7 ACKNOWLEDGEMENT
decided to assign different weights (𝑊 ) to each timestamp while We are grateful to the anonymous reviewers for their constructive
maintaining weight-tied for 𝑉 . Alternatively, the model comprises comments and suggestions. This work was supported in part by
𝑇 layers of GCN and one linear layer. The linear layer serves to the NSF Cybertraining 2320980, SCH 2306331, CAREER 1844403
aggregate static information, while the GCNs handle dynamic in- and the CDC MInD Healthcare group under cooperative agreement
formation. To validate our architecture choice, we conducted a U01-CK000594.
thorough comparison of our model against other configurations:
• Share both: Both 𝑊 and 𝑉 are shared across layers.
• Share 𝑉 : Only 𝑉 is shared across layers.
• Not share: 𝑊 and 𝑉 are not shared.
• W/o loop: as defined in Section 5.4.
The results are presented in Table 7. As observed, the share-
both model exhibits the poorest performance. We believe this is
due to the limited number of free variables available for learning,
which makes the training process challenging. In our approach and
the share-𝑉 method, we achieve very similar results. Our model
utilizes 𝑇 layers of GCN for dynamic information and one linear

Table 7: Evaluating different variants of our model. Present-

ing AUC results on Brain10 datasets

Share both IDGNN Share 𝑉 Not share W/o loop

AUC 94.29 94.87 94.87 94.90 85.62
Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

REFERENCES [28] T Konstantin Rusch, Michael M Bronstein, and Siddhartha Mishra. 2023. A survey
[1] Shaojie Bai, J Zico Kolter, and Vladlen Koltun. 2019. Deep equilibrium models. on oversmoothing in graph neural networks. arXiv preprint arXiv:2303.10993
Advances in Neural Information Processing Systems 32 (2019). (2023).
[2] Qi Chen, Yifei Wang, Yisen Wang, Jiansheng Yang, and Zhouchen Lin. 2022. [29] Aravind Sankar, Yanhong Wu, Liang Gou, Wei Zhang, and Hao Yang. 2020.
Optimization-induced graph implicit nonlinear diffusion. In International Confer- Dysat: Deep neural representation learning on dynamic graphs via self-attention
ence on Machine Learning. PMLR, 3648–3661. networks. In Proceedings of the 13th international conference on web search and
[3] Benoît Colson, Patrice Marcotte, and Gilles Savard. 2007. An overview of bilevel data mining. 519–527.
optimization. Annals of operations research 153 (2007), 235–256. [30] Youngjoo Seo, Michaël Defferrard, Pierre Vandergheynst, and Xavier Bresson.
[4] Laurent El Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Tsai. 2018. Structured sequence modeling with graph convolutional recurrent net-
2021. Implicit deep learning. SIAM Journal on Mathematics of Data Science 3, 3 works. In Neural Information Processing: 25th International Conference, ICONIP
(2021), 930–958. 2018, Siem Reap, Cambodia, December 13-16, 2018, Proceedings, Part I 25. Springer,
[5] Jianfei Gao and Bruno Ribeiro. 2022. On the equivalence between temporal and 362–373.
static equivariant graph representations. In International Conference on Machine [31] Mengzhu Sun, Xi Zhang, Jiaqi Zheng, and Guixiang Ma. 2022. Ddgcn: Dual
Learning. PMLR, 7052–7076. dynamic graph convolutional networks for rumor detection on social media. In
[6] Zhengyang Geng, Xin-Yu Zhang, Shaojie Bai, Yisen Wang, and Zhouchen Lin. Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 36. 4611–4619.
2021. On training implicit models. Advances in Neural Information Processing [32] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero,
Systems 34 (2021), 24247–24260. Pietro Liò, and Yoshua Bengio. 2018. Graph Attention Networks.
[7] Fangda Gu, Heng Chang, Wenwu Zhu, Somayeh Sojoudi, and Laurent El Ghaoui. arXiv:1710.10903 [[Link]]
2020. Implicit graph neural networks. Advances in Neural Information Processing [33] Zehong Wang, Qi Li, and Donghua Yu. 2022. TPGNN: Learning High-order
Systems 33 (2020), 11984–11995. Information in Dynamic Graphs via Temporal Propagation. arXiv preprint
[8] William L. Hamilton, Rex Ying, and Jure Leskovec. 2018. Inductive Representation arXiv:2210.01171 (2022).
Learning on Large Graphs. arXiv:1706.02216 [[Link]] [34] Dongkuan Xu, Wei Cheng, Dongsheng Luo, Yameng Gu, Xinyu Liu, Jingchao Ni,
[9] Danfeng Hong, Lianru Gao, Jing Yao, Bing Zhang, Antonio Plaza, and Joce- Bo Zong, Haifeng Chen, and Xiang Zhang. 2019. Adaptive Neural Network for
lyn Chanussot. 2020. Graph convolutional networks for hyperspectral image Node Classification in Dynamic Networks. 2019 IEEE International Conference on
classification. IEEE Transactions on Geoscience and Remote Sensing 59, 7 (2020), Data Mining (ICDM) (2019), 1402–1407.
5966–5978. [35] Da Xu, Chuanwei Ruan, Evren Korpeoglu, Sushant Kumar, and Kannan Achan.
[10] Quanqi Hu, Yongjian Zhong, and Tianbao Yang. 2022. Multi-block min-max 2020. Inductive representation learning on temporal graphs. arXiv preprint
bilevel optimization with applications in multi-task deep auc maximization. arXiv arXiv:2002.07962 (2020).
preprint arXiv:2206.00260 (2022). [36] Menglin Yang, Ziqiao Meng, and Irwin King. 2020. Featurenorm: L2 feature nor-
[11] Shima Khoshraftar and Aijun An. 2022. A survey on graph representation malization for dynamic graph embedding. In 2020 IEEE International Conference
learning methods. arXiv preprint arXiv:2204.01855 (2022). on Data Mining (ICDM). IEEE, 731–740.
[12] Thomas N Kipf and Max Welling. 2016. Semi-supervised classification with graph [37] Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and
convolutional networks. arXiv preprint arXiv:1609.02907 (2016). Haifeng Li. 2019. T-gcn: A temporal graph convolutional network for traffic
[13] Fuxian Li, Jie Feng, Huan Yan, Guangyin Jin, Fan Yang, Funing Sun, Depeng Jin, prediction. IEEE transactions on intelligent transportation systems 21, 9 (2019),
and Yong Li. 2023. Dynamic graph convolutional recurrent network for traffic 3848–3858.
prediction: Benchmark and solution. ACM Transactions on Knowledge Discovery
from Data 17, 1 (2023), 1–21.
[14] Mingjie Li, Yifei Wang, Yisen Wang, and Zhouchen Lin. 2022. Unbiased Stochastic
Proximal Solver for Graph Neural Networks with Equilibrium States. In The
Eleventh International Conference on Learning Representations.
[15] Qimai Li, Zhichao Han, and Xiao-Ming Wu. 2018. Deeper insights into graph
convolutional networks for semi-supervised learning. In Proceedings of the AAAI
conference on artificial intelligence, Vol. 32.
[16] Ruoyu Li, Sheng Wang, Feiyun Zhu, and Junzhou Huang. 2018. Adaptive Graph
Convolutional Neural Networks. arXiv:1801.03226 [[Link]]
[17] Jie Liao, Wei Zhou, Fengji Luo, Junhao Wen, Min Gao, Xiuhua Li, and Jun Zeng.
2022. SocialLGN: Light graph convolution network for social recommendation.
Information Sciences 589 (2022), 595–607.
[18] Juncheng Liu, Kenji Kawaguchi, Bryan Hooi, Yiwei Wang, and Xiaokui Xiao.
2021. Eignn: Efficient infinite-depth graph neural networks. Advances in Neural
Information Processing Systems 34 (2021), 18762–18773.
[19] Yao Ma, Ziyi Guo, Zhaocun Ren, Jiliang Tang, and Dawei Yin. 2020. Stream-
ing graph neural networks. In Proceedings of the 43rd international ACM SIGIR
conference on research and development in information retrieval. 719–728.
[20] George Panagopoulos, Giannis Nikolentzos, and Michalis Vazirgiannis. 2021.
Transfer Graph Neural Networks for Pandemic Forecasting. In Proceedings of the
35th AAAI Conference on Artificial Intelligence.
[21] Aldo Pareja, Giacomo Domeniconi, Jie Chen, Tengfei Ma, Toyotaro Suzumura,
Hiroki Kanezashi, Tim Kaler, Tao Schardl, and Charles Leiserson. 2020. Evolvegcn:
Evolving graph convolutional networks for dynamic graphs. In Proceedings of
the AAAI conference on artificial intelligence, Vol. 34. 5363–5370.
[22] Junyoung Park, Jinhyun Choo, and Jinkyoo Park. 2021. Convergent Graph Solvers.
In International Conference on Learning Representations.
[23] Junhui Park, Gaeun Sung, SeungHyun Lee, SeungHo Kang, and ChunKyun Park.
2022. ACGCN: graph convolutional networks for activity cliff prediction between
matched molecular pairs. Journal of Chemical Information and Modeling 62, 10
(2022), 2341–2351.
[24] Barak A Pearlmutter. 1994. Fast exact multiplication by the Hessian. Neural
computation 6, 1 (1994), 147–160.
[25] Zi-Hao Qiu, Quanqi Hu, Yongjian Zhong, Lijun Zhang, and Tianbao Yang. 2022.
Large-scale stochastic optimization of ndcg surrogates for deep learning with
provable convergence. arXiv preprint arXiv:2202.12183 (2022).
[26] Emanuele Rossi, Ben Chamberlain, Fabrizio Frasca, Davide Eynard, Federico
Monti, and Michael Bronstein. 2020. Temporal graph networks for deep learning
on dynamic graphs. arXiv preprint arXiv:2006.10637 (2020).
[27] Halsey Lawrence Royden and Patrick Fitzpatrick. 1968. Real analysis. Vol. 2.
Macmillan New York.
Efficient and Effective Implicit Dynamic Graph Neural Network Conference’17, July 2017, Washington, DC, USA

A EXPERIMENT   𝑀2 ... 

0
 
𝑀1 ˆ .. ˜
A.1 Experiment Setup Let 𝑀 := ˆ where 𝑀 :=  0 0  . Let 𝑀 :=
 
𝑀 0  .
0 ... 𝑀𝑡 
We follow the procedure from [5] and utilize the provided code 
𝑀ˆ
 
as the code base to compare all the baselines and our method. 0
. Then
Specifically, we first split the dataset into three portions with ratios 0 𝑀1
70%-10%-20% for training, validation, and testing, respectively. Split- 
0 𝐼𝑚
 
0 𝐼𝑚

ting is based on nodes for transductive tasks and time for inductive ∥|𝑀 |∥𝑜𝑝 = ∥|𝑀˜ | ∥ ≤ ∥|𝑀˜ |∥𝑜𝑝 · ∥ ∥
𝐼𝑛 0 𝑜𝑝 𝐼𝑛 0 𝑜𝑝
tasks. We then normalize node attributes and edge attributes with
the 0-1 normalization method. We train on the training portion, = ∥|𝑀˜ |∥𝑜𝑝 = max{∥|𝑀1 |∥𝑜𝑝 , ..., ∥|𝑀𝑡 |∥𝑜𝑝 }
find the best hyperparameters using the validation set, and report This means if all subsystems satisfy the largest eigenvalue con-
the performance on the test set. We also use ROCAUC score to straint, then the coupled equilibrium equation has a fixed point
evaluate classification tasks and mean average percentage error solution by Lemma A.2. □
(MAPE) for regression tasks.
Proof of Theorem 3.2
A.2 Hyperparameter Theorem 3.2. Let 𝜎 be an element-wise non-expansive function.
For detailed baselines’ architecture, please refer to [5]. Notice that, If the coupled equilibrium equations satisfy the well-posedness con-
for all the methods and all task, we fixed the embedding size as dition, namely ∥M 𝑡 ∥𝑜𝑝 ≤ ∥𝑊 𝑡 ∥𝑜𝑝 ∥𝐴𝑡 ∥𝑜𝑝 < 1, ∀𝑡 = 1, ...,𝑇 , then
16, and we searched the learning rates from 0.1, 0.01, and 0.001. there exists rescale coupled equilibrium equations, which satisfy the
For Brain10, we observed that our method converged slowly, then condition ∥𝑊 𝑡 ∥ ∞ ∥𝐴𝑡 ∥𝑜𝑝 < 1, ∀𝑡 = 1, ...,𝑇 , and the solutions of these
we let it run for 1000 epochs. For DBLP5 and Reddit4, we let all two equations are equivalent.
the methods run 500 epochs for 5 times for each learning rate and
Proof. Suppose {𝑊 𝑡 } satisfy ∥𝑊 𝑡 ∥𝑜𝑝 ∥𝐴𝑡 ∥𝑜𝑝 < 1 for all 𝑡, then
report the performance on the test set where the performance of the
the following equation has a unique fixed point.
validation set is the best. For regression datasets, we run 100 epochs
for England-COVID and 10 epochs for PeMS04/PeMS08. The hyper- 𝑧1   0 0 ··· 0 𝑀 1   𝑧 1   vec(𝑉 𝑋 1 ) 
0   𝑧 2   vec(𝑉 𝑋 2 )  ª®
 2 © 𝑀 2    
parameter our model 𝜂 1 ∈ {0.5, 0.7, 0.9, 1}, 𝜂 2 ∈ {0.001, 0.01, 0.1}. 𝑧  ­ 0 · · · 0
0   𝑧  +  vec(𝑉 𝑋 )  ®®
 3 3 3 3
  = 𝜎 ­­  0
𝑧  ­  𝑀 ··· 0
PROOFS  . 
 ..  ­  ..
­ .
..
.
..
.
..
.
..   ..  
.   .  
..
.
 ®®
®
  
Lemma A.1. If ∥.∥𝑜𝑝 is the matrix operator norm on R𝑛×𝑛 and 𝑧𝑇 
«
 0 0 · · · 𝑀 𝑇 0  𝑧𝑇  vec(𝑉 𝑋 𝑇 ) 
     ¬
𝜆PF (𝐴) outputs the largest absolute eigenvalue of 𝐴, then, for any
matrix 𝐴 ∈ R𝑛×𝑛 , and this condition implies ∥M ∥𝑜𝑝 ≤ 1. Based on Theorem 4.3, there
exists a set of diagonal matrices {𝑆 𝑡 } such that
𝜆PF (𝐴) ≤ ∥𝐴∥𝑜𝑝
𝑊ˆ 𝑡 = 𝑆 𝑡 𝑊 𝑡 (𝑆 𝑡 ) −1, 𝑉ˆ = 𝑆 𝑡 𝑉
Proof. Let 𝜆 be the eigenvalue of 𝐴, and let 𝑥 ≠ 0 be the corre- Then the fixed-point of following equations
sponding eigenvector. From 𝐴𝑥 = 𝜆𝑥, we have  𝑧ˆ1   0 0 ··· 0 𝑀ˆ 1   𝑧ˆ1   vec(𝑉ˆ 𝑋 1 ) 
0   𝑧ˆ2   vec(𝑉ˆ 𝑋 2 )  ®®
 2 © ˆ 2    ª
𝐴𝑋 = 𝜆𝑋  𝑧ˆ  ­ 𝑀 0 · · · 0
 3 3
𝑀ˆ 3 0   𝑧ˆ  +  vec(𝑉ˆ 𝑋 )  ®®
3
­
 𝑧ˆ  ···
  = 𝜎 ­ 0 0
­
where each column in 𝑋 is 𝑥. Further, we have  . 
 .. 
­ .
­  .. .
.. . .. .
.. ..   ..   .. ®
  ­ .  .   . ®
®
|𝜆|∥𝑋 ∥𝑜𝑝 = ∥𝜆𝑋 ∥𝑜𝑝 = ∥𝐴𝑋 ∥𝑜𝑝 ≤ ∥𝐴∥𝑜𝑝 ∥𝑋 ∥𝑜𝑝 𝑧ˆ𝑇 
  «
 0 0 ··· 𝑀 ˆ 𝑇 0    vec(𝑉 𝑋 )  ¬
 𝑧ˆ𝑇  ˆ 𝑇

Since ∥𝑋 ∥ > 0, taking the maximal 𝜆 gives the result. □ satisfies the following relation
 𝑧ˆ1   (𝑆 1 ) −1 ··· 0   𝑧ˆ1 
    
 ..   .. .. ..   .. 
Lemma A.2. The equilibrium equation 𝑧 = 𝜎 (𝑀𝑧 +𝑏) has a unique  . = . . .   . 
    
fixed point solution if ∥|𝑀 |∥𝑜𝑝 < 1, where ∥.∥𝑜𝑝 is the operator norm, 𝑧ˆ𝑇   0
   ··· (𝑆𝑇 ) −1  𝑧ˆ𝑇 
 
and 𝜎 (·) is an element-wise non-expansive function. □

Proof. Based on Lemma A.1 and Theorem 4.1 in [7], this lemma Proof of Remark 1
is immediately followed. □ 𝑁 , we can construct
Proof. Given a set of dynamic graphs {G𝑖 }𝑖=1
a single dynamic graph by merging the snapshots that are from the
Proof of Theorem 3.1 same time stamp, then we obtain
Proof. By Lemma. A.2, the well-posedness requirement for For- Ĝ = {[𝐺𝑖1, ..., 𝐺 𝑁
1
], ..., [𝐺𝑖𝑇 , ..., 𝐺𝑇𝑁 ]} (9)
mula. (3) is ∥M ∥𝑜𝑝 | ≤ 1. Since Formula. (3) and (2) are equivalent,
the well-posedness requirement for Formula. (2) is also ∥|M|∥𝑜𝑝 < Let 𝐴ˆ𝑡
denote the adjacency matrix of [𝐺𝑖𝑡 , ..., 𝐺 𝑁
By theorem 𝑡 ].

1. 3.2, we need to ensure ∥𝑊 𝑡 ∥ ∞ ∥𝐴ˆ𝑡 ∥𝑜𝑝 < 1. Since 𝐴ˆ𝑡 contains

Conference’17, July 2017, Washington, DC, USA Yongjian Zhong, Hieu Vu, Tianbao Yang, and Bijaya Adhikari

Layers GCN-GRU T-GCN IDGNN

8 0.6217 0.9934 1.1113
16 0.5204 0.8719 1.0982
32 0.0077 0.7176 1.0019
Table 8: Smoothness of embeddings. (The larger the better)

𝑁 disconnected graphs, ∥𝐴ˆ𝑡 ∥𝑜𝑝 ≤ max𝑖 ∥𝐴𝑖𝑡 ∥𝑜𝑝 , which means

∥𝑊 𝑡 ∥ ∞ ∥𝐴𝑖𝑡 ∥𝑜𝑝 < 1 needed to be satisfied for all 𝑖. Since 𝑡 is ar-
bitrary, the remark holds. □

B EMBEDDING VISUALIZATION
In this section, we explore an interesting aspect of our method (a) Our model
that can provide empirical insights into its ability to mitigate over-
smoothing. We conduct experiments on a synthetic dataset that
bears resemblance to toy datasets.
The dataset comprises 10 snapshots, with each snapshot rep-
resenting a clique of 10 nodes. Each node is associated with 10
attributes. The nodes fall into two distinct classes, but we delib-
erately conceal the label information in the attributes of the first
snapshot. Specifically, the first two dimensions of the attributes
represent the one-hot encoding of the labels, while the remaining
dimensions are set to zero. Additionally, we assign unique IDs to
the nodes in sequential order. Nodes with IDs ranging from 0 to 4
belong to class 0, while those with IDs ranging from 5 to 9 belong
to class 1. To assess the effectiveness of our method, we visually
compare the embedding results with those obtained from TGCN. (b) TGCN
Upon examining the visualizations, we observe that our model’s
embeddings exhibit gradual changes, whereas TGCN’s embeddings Figure 3: The embedding visualization of our method and
remain consistent for nodes belonging to the same class. From a TGCN
node-centric perspective, TGCN’s embeddings seem reasonable.
Nodes of the same class possess identical features and exhibit the
same topological structure. Therefore, it is logical for them to share a
common embedding. However, our embeddings tend to differentiate
each individual node. We believe that this characteristic plays a
role in mitigating the oversmoothing problem within our model.
Furthermore, we conducted an additional experiment to quanti-
tatively evaluate our model’s ability to tackle over-smoothing. This
experiment is conducted on the binary toy dataset: the toy data
we constructed consists of a dynamic graph with a maximum of
64 snapshots (adapt to layers), with each snapshot being a clique
of 10 nodes. Each node has 10 associated attributes. The task is
binary classification where each node’s class information is hidden
in its first time-stamp’s attributes. Attributes of other time stamps
are randomly sampled from the normal distribution. All methods
are trained with a maximum of 2000 epochs and a learning rate Figure 4: Additional result.
of 0.001. Finally, we evaluate their smoothness using the Mean
Average Distance (MAD). Results are summarized in Table 8.
D ADDITIONAL RESULT ON TOY DATA
C HESSIAN VECTOR PRODUCT This synthetic experiment shifts label information from time stamp
𝜕2 𝑓 1 to time stamp 5. This adjustment ensures uniform difficulty in
To compute the product of Hassian and a vector: 𝐻𝑣, and 𝐻 = 𝜕𝑥 2
.
𝜕𝑓
utilizing label information across all models. Implementing this
𝜕( )𝑇 𝑣 change postpones our method towards achieving 100% accuracy
We compute the product by 𝐻𝑣 = 𝜕𝑥 𝜕𝑥 . In this way, we are not
explicitly computing the Hessian. 6 by approximately 50 epochs. However, even with this modification,
the baselines still struggle to fit the data.
6 [Link]