A Fractional Graph Laplacian Approach to Oversmoothing
Authors: Sohir Maskey
Raffaele Paolino
Aras Bacho
Gitta Kutyniok
Venue: NeurIPS 2023 poster (2023)
Keywords: Graph Neural Networks
Graph Neural ODE
Fractional Laplacian
Oversmoothing
TL;DR: We introduce the fractional graph Laplacian neural ODE, which alleviates oversmoothing and is well-suited for both directed and undirected graphs, as well as various levels of homophily.

Abstract: Graph neural networks (GNNs) have shown state-of-the-art performances in various applications. However, GNNs often struggle to capture long-range dependencies in graphs due to oversmoothing. In this paper, we generalize the concept of oversmoothing from undirected to directed graphs. To this aim, we extend the notion of Dirichlet energy by considering a directed symmetrically normalized Laplacian. As vanilla graph convolutional networks are prone to oversmooth, we adopt a neural graph ODE framework. Specifically, we propose fractional graph Laplacian neural ODEs, which describe non-local dynamics. We prove that our approach allows propagating information  between distant nodes while maintaining a low probability of long-distance jumps. Moreover, we show that our method is more flexible with respect to the convergence of the graph’s Dirichlet energy, thereby mitigating oversmoothing. We conduct extensive experiments on synthetic and real-world graphs, both  directed and undirected, demonstrating our method’s versatility across diverse graph homophily levels.  Our
code is available at https://github.com/RPaolino/fLode

A Fractional Graph Laplacian Approach to Oversmoothing
Sohir Maskey"
Raffaele Paolino"
Department of Mathematics,
Department of Mathematics & MCML,
LMU Munich
LMU Munich
maskey @ mmth.lmu.de
paolino @ www.emaalllm.....
Aras Bacho
Gitta Kutyniok
Department Mathematics,
Department of Mathematics & MCML,
LMU Munich
LMU Munich
Abstract
Graph neural networks (GNNs) have shown state-of-the-ar performances in various applications. However, GNNs often struggle to capture long-range dependencies in graphs due to oversmoothing. In this paper, generalize the concept of oversmoothing from undirected to directed graphs. To this aim, we extend the notion of Dirichlet energy by considering a directed symmetrically normalized Laplacian. As vanilla graph convolutional networks are prone to oversmooth, we adopt neural graph ODE framework. Specifically, we propose fractional graph Laplacian neural ODEs, which describe non-local dynamics. We prove that our approach allows propagating information between distant nodes while maintaining a low probability of long-distance jumps. Moreover, we show that our method is more flexible with respect to the convergence of the graph's Dirichlet energy, thereby mitigating oversmoothing. We conduct extensive experiments on synthetic and real-world graphs, both directed and undirected, demonstrating our method's versatility across diverse graph homophily levels. Our code is available on GitHub.
1 Introduction
Graph neural networks (GNNs) (Gori et al., 2005; Scarselli et al., 2009; Bronstein et al., 2017) have emerged as a powerful class of machine learning models capable of effectively learning rep- resentations of structured data. GNNs have demonstrated state-of-the- range of applications, including social network analysis (Monti al., 2019), molecular property prediction (Gilmer et al., 2017), and recommendation systems (J. Wang et al., 2018; Fan et al., 2019). The majority of existing work on GNNs has focused on undirected graphs (Defferrard al., 2016; Kipf al., 2017; Hamilton et al., 2017), where edges have no inherent direction. However, many real-world systems, such as citation networks, transportation systems, and biological pathways, are inherently directed, necessitating the development of methods explicitly tailored to directed graphs.
performance in wide
Despite their success, most existing GNN models struggle capture long-range dependencies, which can be critical for specific tasks, such as node classification and link prediction, and for specific graphs, such as heterophilic graphs. This shortcoming also arises from the problem oversmoothing, where increasing the depth of GNNs results in the node features converging to similar values that only convey information about the node's degree (Oono al., 2019; Cai et al., 2020). Consequently,
"Equal contribution.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

scaling the depth of GNNs is not sufficient to broaden receptive fields, and other approaches are necessary to address this limitation. While these issues have been extensively studied undirected graphs (Q. Li et al., 2018; G. Li et al., 2019; Luan, M. Zhao, et al., 2019; D. Chen et al., 2020; Rusch al., 2022), their implications for directed graphs remain largely unexplored. Investigating these challenges and developing effective solutions is crucial for applying GNNs to real- scenarios.
Over-smoothing has been shown be intimately related to the graph's Dirichlet energy, defined as
N
X; Vd
Xj dd
(
=
f,f=1
12
where A (at))) represents the adjacency matrix of the underlying graph, x E IR N denotes the node features, and d; E the degree of node i. Intuitively, the Dirichlet energy measures the smoothness of nodes' features. Therefore, GNN that minimizes the Dirichlet energy is expected to perform well on homophilic graphs, where similar nodes are likely to be connected. Conversely, a GNN that ensures high Dirichlet energy should lead to better performances heterophilic graphs, for which the nodes' features are less smooth.
This paper aims to bridge the gap in understanding oversmoothing for directed graphs. To this aim, we generalize the concept of Dirichlet energy, providing rigorous foundation for analyzing oversmoothing. Specifically, we consider the directed symmetrically normalized Laplacian, which accommodates directed graph structures and recovers the usual definition in the undirected case. Even though the directed symmetrically normalized Laplacian has been already used (Zou et al., 2022), its theoretical properties remain widely unexplored.
However, a vanilla graph convolutional network (GCN) (Kipf et al., 2017) implementing this directed Laplacian alone is not able to prevent oversmoothing. For this reason, we adopt a graph neural ODE framework, which has been shown to effectively alleviate oversmoothing in undirected graphs (Bodnar et al., 2022; Rusch al., 2022; Di Giovanni al., 2023).
1.1 Graph Neural ODEs
The concept of neural ODE was introduced by Haber et al. (2018) and R. T. Q. Chen et al. (2018), who first interpreted the layers in neural networks the time variable in ODEs. Building on this foundation, Poli et al. (2021), Chamberlain et al. (2021), and Eliasof et al. (2021) extended the connection the realm of GNNs, resulting in the development of graph neural ODEs. In this context, each node of the underlying graph is described by a state variable x; (t) E RR representing the node i time f. We can define the dynamics of x(f) via the node-wise ODE
x' (t) = (x()))) t E [0,T], subject to the initial condition x0) parametrized by the learnable parameters w.
Xo E R NxK where the function fw RNK
R NxK is
The graph neural ODE can be seen as a continuous learnable architecture on the underlying graph, which computes the final node representation x(T) from the input nodes' features Xg- Typical choices for fw include attention-basee networks (GATs) (Velickovic et al., 2018), or convolutional- like functions (Di Giovanni al., 2023) that generalize GCNs (Kipf et al., 2017).
functions (Chamberlain et al., 2021), which generalize graph attention
How can we choose the learnable function fw accommodate both directed and undirected graphs, as well as different levels of homophily? We address this question in the following subsection.
1.2 Fractional Laplacians
The continuous fractional Laplacian, denoted by (-A)a a for 0, is used to model non-local interactions. For instance, the fractional heat equation 0,u + (-)) 0 provides a flexible and accurate framework for modeling anomalous diffusion processes. Similarly, the fractional diffusion- reaction, quasi-geostroph geostrophi,,, Sahn-Hilliard, porous medium, Schrödinger, and ultrasound equations more sophisticated models represent complex anomalous systems (Pozrikidis, 2018).
Similarly to the continuous case, the fractional graph I Laplacian (FGL) (Benzi et al., 2020) models non-local network dynamics. In general, the FGL does not inherit the sparsity of the underlying

---

graph, allowing a random walker to leap rather than walk solely between adjacent nodes. Hence, the FGL able build long-range connections, making it well-suited for heterophilic graphs.
1.3 Main Contributions
We present novel approach to the fractional graph Laplacian by defining it in the singular value domain, instead of the frequency domain (Benzi al., 2020). This formulation bypasses the need for computing the Jordan decomposition of the graph Laplacian, which lacks reliable numerical methods. We show that our version of the FGL can still capture long-range dependencies, and prove that its entries remain reasonably bounded.
We then propose two FGL-based neural ODEs: the fractional heat equation and the fractional Schrödinger equation. Importantly, we demonstrate that solutions to these FGL-based neural ODEs offer increased flexibility in terms of the convergence of the Dirichlet energy. Notably, the exponent of the fractional graph Laplacian becomes a learnable parameter, allowing our network to adaptively determine the optimal exponent for the given task and graph. We show that this can effectively alleviate oversmoothing in undirected and directed graphs.
To validate the effectiveness of our approach, we conduct extensive experiments synthetic and real-world graphs, with a specific focus on supervised node classification. Our experimental results indicate the advantages offered by fractional graph Laplacians, particularly in non-homophilic and directed graphs.
Preliminaries
We denote graph = E), where V is the set of nodes, & is the set of edges, and N |>]| is the number of nodes. The adjacency matrix == {as.j} encodes the edge information, with 0.j = if there is edge directed from node i, and 0 otherwise. The in- and out-degree matrices then defined as Din = diag(A11
out = diag( A 1), respectively. The node feature matrix x E IR N contains for every node its feature in IR
Given any matrix M E Caxn, we denote its spectrum by A(M) = {A,(M))² in ascending order w.r.t. to the real part, i.e., RA1 (M) < (M) ... (M). Furthermore, denote by |MM]] and || |M] the Frobenius and spectral norm of M, respectively. Lastly, we denote by In the identity matrix, where omit the dimension when it is clear from the context.
Homophily and Heterophily Given a graph G E) with labels y = {yi hey. the homophily of the graph indicates whether connected nodes are likely to have the same labels; formally,
1 N
o. = 11 -,
(G)
KiE (1, (1,..., NJ Oi.j 1}
where the numerator represents the number of neighbors of node E V that have the same label y: (Pei al., 2019). We say that G is homophilic if H (g)) xe 1 and heterophilic if H (G) Rs 0.
3 Dirichlet Energy and Laplacian for (Directed) Graphs
In this section, we introduce the concept of Dirichlet energy and demonstrate its relationship to directed Laplacian, thereby generalizing Definition 3.1. The Dirichlet energy is defined on the node features NRK of graph G as
results for undirected graphs
X;
8 (x))
Ø..j
i,j=1
The Dirichlet energy measures how much the features change over the nodes of G, by quantifying the disparity between the normalized outflow of information from node and the normalized inflow of information node i.
Definition 3.2. We define the symmetrically normalized adjacency (SNA) as L = D AD 1/2
AD
our
3

---

Chameleon
Squirrel
3A
3
RA
RA
Figure I: Spectrum À (L) common directed real-world graphs. The Berron-Frobente cigenvalue App se 0.94 for Chameleon, and Xpp 0.89 for Squirrel.
Figure 2: Examples of non-weakly balanced (left), weakly balanced (center), and balanced (right) directed graphs. The Perron-Frobenius eigenvalue of the left graph is Apr 0.97 o 1, while for the middle and right graphs Xpp 1.
Note that L symmetric if and only is undirected; the term "symmetrically" refers to the both-sided normalization rather than the specific property of the matrix itself.
is well-known that the SNA's spectrum of connected undirected graph lies within [-1, 1] (Chung, 1997). We extend this result to directed graphs, which generally exhibit complex-valued spectra. Proposition 3.3. Let be a directed graph with SNA For every À E X(L), it holds | S
Proposition 3.3 provides upper bound for the largest eigenvalue of any directed graph, irrespective of its size. However, many other spectral properties do not carry over easily from the undirected to the directed case. For example, the SNA may possess one-eigenvalue, even if the graph is strongly connected (see, e.g., Figures to 2). The one-eigenvalue is of particular interest since its cigenvector corresponds to zero Dirichlet energy ((v) O. Therefore, studying when 1 E A(L) is crucial to understanding the behavior of the Dirichlet energy. We fully characterize the set of graphs for which XL)) this is the scope of the following definition.
Definition 3.4. A graph 9 (V,E) is said to be balanced if d. det for all E (1,.... N), and weakly balanced if there exists k E RRN such that k + 0 and
hi Va
Σ
aij
0, Vi (1,...
d
It is straightforward that balanced graph is weakly balanced since Hence, all undirected graphs are also weakly balanced. However, as shown in Figure 2, the set of balanced graphs is proper subset of the of weakly balanced graphs.
choose k; / d
Proposition 3.5. Let G be directed graph with SNA L. Then, E X(L) if and only the graph is weakly balanced. Suppose the graph strongly connected, then - X(L) if and only if the graph weakly balanced with even period.
Proposition 3.5 generalizes a well-known result for undirected graphs: - E A(L) if and only if the graph is bipartite, i.e., has even period. The next result shows that the Dirichlet energy defined in (1) and the SNA are closely connected.
Proposition 3.6. For every x E CN holds ((x) R (trace (I- Moreover, there exists # 0 such that & (x) = 0 if and only if the graph is weakly balanced.
Proposition 3.6 generalizes the well-known result from the undirected (see, e.g., Cai et al., 2020, Definition 3. to the directed case. This result an important tool for analyzing the evolution the Dirichlet energy in graph neural networks.

---

20
Q= 2
Q 2
Q=2-
0.5
00.5
(a) Synthetic cycle graph. The values of the fractional Laplacian also be negative. Cora
Chameleon
106
N. of virtual edges
Q
10c
0.75
10a
0.5
10a
0.25
0
10
10
5
10
15
10
15
20
Distance in the original graph
Distance in the original graph
(b) Real world graphs. count number virtual edges built by the fractional Laplacian based distance d(i, in the original graph. The number of virtual edges increases as a decreases.
Figure 3: Visual representation of long-range edges built by the fractional Laplacian.
4 Fractional Graph Laplacians
We introduce the fractional graph Laplacian through the singular value decomposition (SVD). This approach has two key advantages over the traditional definition (Pozrikidis, 2018; Benzi et al., 2020) in the spectral domain. First, it allows defining the fractional Laplacian based on any choice of graph Laplacian, including those with negative or complex spectrum such as the SNA. Secondly, the SVD computationally more efficient and numerically more stable than the Jordan decomposition, which would be necessary if the fractional Laplacian was defined in the spectral domain.
Consider directed graph with SNA L and its SVD L UWVH, where U, V E CN*N are unitary matrices and E R. X is a diagonal matrix. Given a E R, we define the a-fractional graph Laplacian (a-FGL in short)
L = V
In undirected graphs, the a-FGL preserves the sign of the eigenvalues À of while modifying their magnitudes, i.e., À + sign(A) "
The a-FGL is generally less sparse than the original SNA, connects nodes that are not adjacent in the underlying graph. The next theorem proves that the weight of such "virtual" edges is bounded.
Theorem 4.1. Let G be directed graph with SNA L. For a > if the distance d(i,j) between nodes and j is at least 2, then
n² 2 11
L
(L (L)) 1
2 (d
We provide proof of Theorem 4.1 in Appendix C. In Figure 3a, visually represent the cycle graph with eight nodes and the corresponding a-FGL entries. We also refer to Figure 3b, where we
We employ the term "fractional graph Laplacian" instead of "fractional symmetrically normalized adjacency" to highlight that, in the singular value domain, one can construct fractional powers of any matrix representations of the underlying graph (commonly referred to "graph Laplacians"). This fact is not generally true in the eigenvalue domain because eigenvalues can be negative or complex numbers.
"In fact, this property can be generalized to all directed graphs with a normal SNA matrix, including, among others, directed cycle graphs. See Lemma C.I for more details.

---

depict the distribution of a-FGL entries for the real-world graphs Cora (undirected) and Chameleon (directed) with respect to the distance in the original graph. Our empirical findings align with our theoretical results presented in Theorem 4.1.
5 Fractional Graph Laplacian Neural ODE
This section explores two fractional Laplacian-barr based graph neural ODEs. First, we consider the fractional heat equation,
x()) = (t), x(0) = X0,
(2)
where Xo E N is the initial condition, x(t) E RN* for > and a E R. We assume that the channel mixing matrix W E R is symmetric matrix. Second, we consider the fractional Schrödinger equation,
xt) i x(0) Xo
where Xo, x(t E CN and W CK*K is unitary diagonalizable. Both (2) and (3) can be analyt- ically solved. For instance, the solution of (2) is given by vec(x)(t) = exp(-t W & L*)) ve)( where @ denotes the Kronecker product and vec(-) represents the vectorization operation. How- ever, calculating the exact solution is computationally infeasible since the memory required to store W & L alone grows (NK) Therefore, rely numerical schemes to solve (2) and (3).
In the remainder of this section, we analyze the Dirichlet energy for solutions to (2) and (3). We begin with the definition of oversmoothing.
Definition 5.1. Neural ODE-based GNNs are said to oversmooth if the normalized Dirichlet energy decays exponentially fast. That for any initial value Xo, the solution x(t) satisfies for every t
x(1) |x(0]]] [x(t)ll2
min A(I L) S exp (-Ct) C> C> 0.
Definition 5.1 captures the actual smoothness of features by considering the normalized Dirichlet energy, which mitigates the impact of feature amplitude (Cai et al., 2020; Giovanni et al., 2023). Additionally, Proposition 3.6 shows that the normalized Dirichlet energy is intimately related to the numerical range of I- L of the underlying graph. This shows that the Dirichlet energy and cigenvalues (or frequencies) of the SNA are intertwined, and one can equivalently talk about Dirichlet energy or frequencies (see also Lemma D.2). In particular, holds that
x(t)
[ [-LI Lll 2
0 &
|x(t)
As seen in Section 3, the minimal possible value attained by the normalized Dirichlet energy is often strictly greater than 0 for directed graphs. This indicates that GNNs on general directed graphs inherently cannot oversmooth to the same extent as undirected. However, we prove that a vanilla GCN implementing the directed SNA oversmooths with respect to Definition 5.1, see Appendix E.3.
5.1 Frequency Analysis for Graphs with Normal SNA
This subsection focuses on the frequency analysis of FGL-based Neural ODEs for undirected graphs Most classical GNNs (Kipf et al., 2017; Velickovic et al., 2018) and also graph neural ODEs (Chamberlain et al., 2021; Eliasof et al., 2021) have been shown to oversmooth. Di Giovanni et al. (2023) proved that the normalized Dirichlet energy for GNNs based on (2) with = 1 can not only converge to its minimal value but also to its maximal possible value. A GNN exhibiting this property is then termed FiighestFrreuuncyy -Dominnnn (HFD).
However, in real-world scenarios, most graphs are not purely homophilic nor purely heterophilic but fall somewhere in between. Intuitively, this suggests that mid-range frequencies might be more suitable. To illustrate this intuition, consider the cycle graph as an example. If we have homophily of 1, low frequencies are optimal; with homophily equal to 0, high frequencies are optimal. Interestingly, for homophily of 1/2, the mid-range frequency is optimal, even though the eigendecomposition is label-independent. More information on this example can be found in Figure 4 and Appendix F. Based on this observation, we propose the following definition to generalize the concept of HFD, accommodating not only the lowest or highest frequency but all possible frequencies.

---

-=-1
X -=122/2
=
== v2/2
A ==1
0
-22/2
H
Hl 1/2
H
Figure 4: Eigendecomposition of L for the cycle graph Cg (see Appendix F). The first two rows show the eigenvectors corresponding to the eigenvalues A. The last row shows how the (lahel-unaware) eigendecomposittinn can be used to study homophily, whose definition requires the labels.
Definition 5.2. Let O. Neural ODE-based GNNs initialized at Xy are X-Frequency-Domme Dommnnt (A-FD) if the solution x(t) satisfies
x(t)
[-+00 A
(x(2)ll2
Suppose À the smallest or the largest eigenvalue with respect to the real part. In that case, we call it Lowest- Frevesttrrqueeny--DDmmoee Dominant (LFD) gighest-irequenc--Donr -oonnee----Dmntttttoee (HFD), respectively.
In the following theorem, we show that (2) and (3) are not limited being LFD or HFD, but can also be mid-frqquency dominant.
Theorem 5.3. Let G be an undirected graph with SNA L. Consider the initial value problem in (2) with W R and a E R. Then, for almost all initial values Xg RNNK the following holds.
(a 0) The solution to (2) either HFD or LFD.
(a 0) Let (L) and (L) be the smallest positive and negative non-zero eigenvalue of L. respectively. The solution to (2) is either (1 -1- X4 (L))-FD or (1-X (L))-D
Furthermore, the previous results also hold for solutions to the Schrödinger equation (3) if W E CK*K has at least one eigenvalue with non-zero imaginary part.
Theorem 5.3 (a > 0) generalizes the result by Di Giovanni et al. (2023) for positive values a. The convergence speed Thcorem 5.3 (a > 0) depends on the choice of E R. By selecting a variable (c.g., as a learnable parameter), framework capable of adapting the convergence speed the Dirichlet energy. A slower or more adjustable convergence speed facilitates broader frequency exploration converges to its maximal or minimal value. Consequently, the frequency component contributions (for finite time, i.e., in practice) are better balanced, which is advantageous for graphs with different homophily levels. Theorem 5.3 (a < 0) shows that solutions of the fractional neural ODEs in (2) and (3) are not limited to be LFD HFD. To demonstrate this and the other results of Theorem 5.3, solve (2) using an explicit Fuler scheme for different choices of a and W on the Cora and Chameleon graphs. The resulting evolution of the Dirichlet energy with respect to time is illustrated in Figure 5. Finally, we refer to Theorem D.5 in Appendix D.I for the full statement and proof of Theorem 5.3.
to arbitrary
establish flexible learning
gradually
Remark 5.4. Theorem 5.3 stated for the analytical solutions of (2) and (3), respectively. As noted in Section 5, calculating the analytical solution is infeasible in practice. However, we show in Appendices D.. to 1).3 that approximations of the solution of (2) and (3) via explicit Euler schemes satisfy the same Dirichlet energy convergence properties if the step size is sufficiently small.
Remark 5.5. Theorem 5.3 can be generalized to all directed graphs with normal SNA, Le., satisfying the condition LL L L. For the complete statement, see Appendix D.I.
7

---

0.5
Dirichlet energy
a
0.75
0.75
0.25
0.5
0.5
0.25
0.25
0
0.1
100
10
10a
100
10a
10
Time
Time
(a) Cora (undirected).
0.5
Dirichlet energy
a
0.75
0.75
0.25
0.5
0.5
0.25
0
0.25
0.1
100
10a
10a
100
10a
10
Time
Time
(b) chamelcon (directed).
Figure 5: Convergence of Dirichlet energy for the solution of equation (2) using explicit Euler scheme with a step size of h = 10-1 We consider different a-FGL. in (2) and choose W as a random diagonal matrix. In the left plot, W has only negative spectrum, while in the right plot, W has only a positive spectrum. The black horizontal line represents the theoretical limit based on Theorem 5.3.
5.2 Frequency Dominance for Directed Graphs
Section 5.1 analyzes the Dirichlet energy in graphs with normal SNA. However, the situation becomes significantly more complex when considering generic directed graphs. In our experiments (see Figure 5), we observe that the solution to (2) and (3) does not necessarily lead to oversmoothing. On the contrary, the solution can be controlled to exhibit either LFD or HFD for a >> 0, and mid- frequency-dominane for a < as proven for undirected graphs in Theorem 5.3, We present an initial theoretical result for directed graphs, specifically in the case of 1.
Theorem 5.6. Let G be directed graph with SNA L. Consider the initial value problem in (2) with diagonal channel mixing matrix W R x and a 1. Suppose X (L) unique. For almost all initial values E RRN X the solution to (2) is either HFD or LFD.
The proof of Theorem 5.6 is given in Appendix E.1. Finally, we refer to Appendix E.2 for the analogous statement and proof when the solution of (2) is approximated via explicit Euler scheme.
Numerical Experiments
This section evaluates the fractional Laplacian ODEs in node classification by approximating (2) and (3) with an explicit Euler scheme. This leads to the following update rules
X x, 1 L W. X L
+ h
for the heat and Schrödinger equation. respectively. In both cases, W, a and h are learnable parameters, is the layer, and Xi is the imitial nodes' feature matrix. In accordance with the results in Section 5, we select W as diagonal matrix. The initial features Xo (4) are encoded through MI.,, and the output is decoded using second MLP. We refer to the resulting model as FLODE (fractional Laplacian ODE). In Appendix A, we present details on the baseline models, the training setup, and the exact hyperparameters.

---

Table 1: Test accuracy (Film. Squirrel, Chameleon, Citeseer) and test AUROC (Minesweeper. Tolokers, Questions) on node classification, top three models. The thorough comparison is reported in Table 4, Appendix A: FLODE consistently outperforms the baseline models GCN and GRAFF, and it achieves results comparable to state-of-the-art.
(a) Undirected graphs.
(b) citeroppill-sppeciffe graphs.
Squirrel
Chameleon
Citeseer
Minesweeper
Tolokers
Questions
FLODE
FLODE 73.60 + 1.55
FLODE 78.07 + 1.62
GeT 93.91 + 0.35
FLODE 84.17 + 0.58
FSGNN 78.86 + 0.92
1"
64.23 + 1.84
2nd
GREAD
GREAD 71.38 + 1.30
Geom-GCN 78.02 + 1.15
GraphSAGE 93.51 + 0.57
SAT 83.78 + 0.43
FLODE 78.39 + 1.22
59.22 + 1.44
3rd
GRAFFN.. 59.01 + 1.31
GRAFFNL 71.38 + 1.47 77.60 + 1.81
GREAD
FLODE 92.43 + 0.51
GAT 83.70 + 0.47
GT-sep 78.05 + 0.93
(c) Directed graphs.
Film
Squirrel
Chamelcon
FLODE
HLP 71.17 + 1.83
FSGNN 78.14 + 1.25
37.41 = 1.06
2nd
GRAFF
FLODE 74.03 + 1.58
FLODE 77.98 + 1.05
37.11 = 1.08
3
ACM
ESGNN 73.48 + 2.13
HLP 77.48 + 1.50
36.89 1.18
Ablation Study. In Appendix A.3, we investigate the influence of each component (learnable exponent, ODE framework, directionality via the SNA) on the performance of FLODE. The adjustable fractional power in the FGL is a crucial component of FLODE, as it alone outperforms the model employing the ODE framework with a fixed 1. Further, Appendix A.3 includes ablation studies that demonstrate FLODE's capability to efficiently scale to large depths, as depicted in Figure 8.
Real-Worle Graphs. We report results on 6 undirected datasets consisting of both homophilic graphs, i.e., Cora (McCallum et al., 2000), Citeseer (Sen et al., 2008) and Pubmed (Namata et al., 2012), and heterophilic graphs, i.e., Film (Tang et al., 2009), Squirrel and Chameleon (Rozemberezki et al., 2021). We evaluate our method on the directed and undirected versions of Squirrel, Film, and Chameleon. In all datasets, we use the standard 10 splits from (Pei et al., 2019). The choice of the bascline models and their results are taken from (Di Giovanni et al., 2023). Further, we test our method on heterophily-specinn Questions (Platonov et al., 2023). The splits, baseline models, and results are taken from (Platonov et al., 2023). The top three models are shown in Table 1, and the thorough comparison is reported in Table 4. Due to memory limitations, compute only 30% of singular values for Pubmed, Roman- Empane, and Questions, which serve as the best low-rank approximation of the original SNA.
graph datasets, i.e., Roman-empire, Minesweeper, Tolokers, and
Synthetic Directed Graph. We consider the directed stochastic block model (DSBM) datasets (Zhang et al., 2021). The DSBM divides nodes into 5 clusters and assigns probabilities for interactions between vertices. It considers two sets of probabilities: {as.s} for undirected edge creation and { for assigning edge directions, i j E {1,...5), 5). The objective is to classify vertices based on their clusters. In the first experiment, æi.j = a varies. altering neighborhood information's 's importance. In the second experiment, Br,j = 3 varies, changing directional information. The results are shown in Figure and Table The splits, baseline models, and results are taken from (Zhang et al., 2021).
The experiments showcase the flexibility of FLODE, as it can accommodate various types of graphs, both directed and undirected, as well as a broad range of homophily levels. While other methods, such as MagNet (Zhang et al., 2021), perform similarly to our approach, they face limitations when applied to certain graph configurations. For instance, when applied to undirected graphs, MagNet reduces to ChebNet, making it unsuitable for heterophilic graphs. Similarly, GRAFF
Results.

---

100
MagNet
DiGraphIB
DiGraph
Accuracy
DGCN
50
GAT-D
GIN-D
GraphSAGE-D
APPNP-D
GCN-D
0.1
0.08
0.05
0.1
0.2
0.3
0.4
ChebNet
o=
8
FLODE
Figure 6: Experiments directed stochastic block model. Unlike other models, FLODEss perfor- mances do not deteriorate as much when changing the inter-cluster edge density a²
Chamelcon
Squirrel
Explained Variance
Explained Variance
10
0.8
Test accuracy
109
0.75
Test accuracy
0.65
0.7
10
10
0.55
0.6
0.45
10
0.5
10
0.35
100
10
102
103
10
10a
102 103
N. singular values
N. of singular values
Figure 7: Efffect of truncated SVD on test accuracy (orange) for standard directed real-world graphs. The explained variance, defined as Et= <?/>;"- of. measures the variability the first & singular values explain. For chameleon, the accuracy stabilized after 570 (25%) singular values, corresponding to an explained variance of 0.998. For squirrel, after 1600 (31 %) singular values, which correspond an explained variance 0.999, the improvement in test accuracy is only marginal.
(Di Giovanni et al., 2023) performs well on undirected graphs but falls short on directed graphs. We note that oftentimes FLODE learns non-trivial exponent 1, highlighting the advantages of FGLebased GNNs (see, e.g., Table 5). Furthermore, as shown in Table and Appendix A.3, our empirical results align closely with the theoretical results in Section
7 Conclusion
In this work, we introduce the concepts of Dirichlet energy and oversmoothing for directed graphs and demonstrate their relation with the SNA. Building upon this foundation, we define fractional graph Laplacians in the singular value domain, resulting in matrices capable of capturing long-range dependencies. To address oversmoothing directed graphs, we propose fractional Laplacian- based graph ODEs, which are provably not limited to LFD behavior. We finally show the flexibility of our method to accommodate various graph structures and homophily levels in node-level tasks.
Limitations and Future Work. The computational cost of the SVD grows cubically in N, while storage the singular vectors grows quadratically in N. Both costs can significantly reduced by computing only k & N singular values via truncated SVD (Figure 7), giving the best k-rank approximation of the SNA. Moreover, the SVD can be computed offline as a preprocessing step.
The frequency analysis of a-FGL neural ODEs in directed graphs is an exciting future direction. It would also be worthwhile to investigate the impact of choosing a of the convergence speed of the Dirichlet energy. Controlling the speed could facilitate the convergence of the Dirichlet energy to an optimal value, which has been shown exist synthetic settings (Keriven, 2022; X. Wu et al., 2022). Another interesting future direction would be to analyze the dynamics when approximating the solution the FGL neural ODEs using alternative numerical solvers, such as adjoint methods.
10

---

Acknowledgments
M. acknowledges partial support by the NSF-Simons Research Collaboration on the Mathematical and Scientific Foundations of Deep Learning (MoDL) (NSF DMS 2031985) and DFG SPP 1798, KU 1446/27-2.
G. K. acknowledges partial support by the DAAD programme Konrad Zuse Schools of Excellence in Artificial Intelligence, sponsored by the Federal Ministry of Education and Research. G. Kutyniok also acknowledges support from the Munich Center for Machine Learning (MCML) as well as the German Research Foundation under Grants DFG-SPP-2298, KU 1446/31-1 and KU 1446/32-1 and under Grant DFG-SFB/TR 109 and Project CO9.
References
Benzi, M., Bertaccini, D., Durastante, F., and Simunec, I. (2020). "Non-Local Network Dynamics via Fractional Graph Laplacians". In: Journal of Complex Networks 8.3.
Bo, D., Wang, X., Shi, C., and Shen, H. (2021). "Beyond Low-frequency Information in Graph Convolutional Networks". In: Proceedings of the AAAI Conference on Artificial Intelligence 35.5, pp. 3950-3957.
Bodnar, C., Di Giovanni, F., Chamberlain, B., Lió, P., and Bronstein, M. (2022). "Neural Sheaf Diffusion: A Topological Perspective on Heterophily and Oversmoothing in GNNs". In: Advances in Neural Information Processing Systems 35, pp. 18527-18541.
Bronstein, M. M., Bruna, J., LeCun, Y.. Szlam, A., and Vandergheynst, P. (2017). "Geometric Deep Learning: Going beyond Euclidean Data". In: IEEE Signal Processing Magazine 34.4, pp. 18-42. Cai, C. and Wang. Y. (2020). A Note on Over-Smoothing for Graph Neural Networks. arxiv: 2006. 13318 [cs, stat].
Chamberlain, B., Rowbottom, J., Gorinova, M. I., Bronstein, M., Webb, S., and Rossi, E. (2021).
"GRAND: Graph Neural Diffusion". In: Proceedings of the 38th International Conference on Machine Learning. PMLR, pp. 1407-1418.
Chen, D., Lin, Y., Li, W., Li, P., Zhou, J., and Sun, X. (2020). " "Measuring and Relieving the Over- Smoothing Problem for Graph Neural Networks from the Topological View". In: Proceedings of the AAAI Conference on Artificial Intelligence 34.04, pp. 1438-3445
Chen, M., Wei, Z., Huang, Z., Ding, B., and Li, Y. (2020). "Simple and Deep Graph Convolutional Networks". In: Proceedings of the 37th International Conference on Machine Learning. PMLR,
pp. 1725-1735.
Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. (2018). "Neural Ordinary Differential Equations". In: Advances in Neural Information Processing Systems. Vol. 31. Curran
Associates, Inc.
Chien, E., Peng, J., Li, P., and Milenkovic, O. (2021). " "Adaptive Universal Generalized PageRank Graph Neural Network". In: International Conference on Learning Representations.
Choi, J., Hong, S., Park, N., and Cho, S.-B. (2023). "GREAD: Graph Neural Reaction-Diffusion
Networks". In: Proceedings of the 40th International Conference on Machine Learning. Vol. 202. ICML'23, Honolulu, Hawaii, USA: JMLR.org, pp. 5722-5747.
Chung. F. R. K. (1997). Spectral Graph Theory. Regional Conference Series in Mathematics no. 92. Providence, R.I: Published for the Conference Board of the mathematical sciences by the American Mathematical Society.
Defferrard, M., Bresson, X., and Vandergheynst, (2016). "Convolutional Neural Networks on Graphs with Fast Localized Spectral Filtering". In: Advances in Neural Information Processing Systems. Vol. 29, Curran Associates, Inc.
Di Giovanni, F., Rowbottom, J., Chamberlain, B. P., Markovich, T., and Bronstein, M. M. (2023). "Understanding convolution on graphs via energies". In: Transactions on Machine Learning Research.
Du, L., Shi, X., Fu, Ma, X., Liu, H., Han, S., and Zhang, D. (2022). "GBK-GNN: Gated Bi-Kernel Graph Neural Networks for Modeling Both Homophily and Heterophily". . In: Proceedings of the ACM Web Conference 2022. WWW '22. New York, NY, USA: Association for Computing Machinery, pp. 1550-1558.
Eliasof, M., Haber, E., and Treister, E. (2021). "PDE-GCN: Novel Architectures for Graph Neural Networks Motivated by Partial Differential Equations". In: Advances in Neural Information Processing Systems. Vol. 34. Curran Associates, Inc., pp. 3836-3849.
11

---

Fan. W., Ma, Li, Q., He, Y., Zhao, E.,. Tang. J., and Yin, D. (2019). "Graph Neural Networks for Social Recommendation" In: The World Wide Web Conference. San Francisco CA USA: ACM, pp. 417-426.
Fey, M. and Lenssen, J. F. (2019). "Fast Graph Representation Learning with PyTorch Geometric".
In: ICLR 2019 Workshop on Representation Learning on Graphs and Manifolds. Gasteiger, J., Bojchevski, A., and Günnemann, S. (2018). "Predict Then Propagate: Graph Neural
Networks Meet Personalized PageRank". In: International Conference on Learning Representa-
tions.
Gilmer, J., Schoenholz, S. S., Riley, P. F., Vinyals, O., and Dahl, G. E. (2017). "Neural Message Passing for Quantum Chemistry". In: Proceedings of the 34th International Conference on Machine Learning. PMLR. pp. 1263-1272
Gori, M., Monfardini, G., and Scarselli, F. (2005). "A New Model for Learning in Graph Domains". In: P Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005. Vol. 2, 729-734 vol. 2.
Haber, E. and Ruthotto, L. (2018). "Stable Architectures for Deep Neural Networks". In: Inverse Problems 34.1, p. 014004, eprint: 1705.03341 (cs, math).
Hamilton, W., Ying, Z., and Leskovec, J. (2017). "Inductive Representation Learning on Large Graphs". In: Advances in Neural Information Processing Systems. Vol. 30. Curran Associates, Inc. He, K., Zhang. X., Ren, S., and Sun, (2016). "Deep Residual Learning for Image Recognition", In: Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 770-778. Keriven, N. (2022). "Not Too Little, Not Too Much: A Theoretical Analysis of Graph (over)Smoothh ing". In: Advances in Neural Information Processing Systems.
Kingma, D. P. and Ba, (2015). Adam: A Method for Stochastic Optimization. Ed. by Y. Bengio and Y. LeCun.
Kipf. T. N. and Welling, M. (2017). "SeemiSSuerriiie Classification with Graph Convolutional Networks". In: International Conference Learning Representations.
Li, G., Muller, M., Thabet, A., and Ghanem, B. (2019). "DeepGCNs: Can GCNs Go As Deep As CNNs?" In: Proceedings of the IEEE/CVF International Conference on Computer Vision, pp. 90699--9999 9267-9276.
Li, Q., Han, Z., and Wu, X.m (2018). "Deeper Insights Into Graph Convolutional Networks for Scmi-Supervised Learning". In: Proceedings of the AAAI Conference on Artificial Intelligence
32.1.
Li, X., Zhu, R., Cheng, Y., Shan, C., Luo, S., Li, D.,, and Qian, W. (2022). "Finding Global Homophily in Graph Neural Networks When Meeting Heterophily". In: Proceedings of the 39th International Conference on Machine Learning. PMLR, pp. 19242-3256.
Lingam, V., Ragesh, R., Iyer, A., and Sellamanickam, S. (2021). Simple Truncated SVD Based Model for Node Classification on Heterophilic Graphs. eprint: 2106 12807 (cs). Luan, S., Hua, C., Lu, Q.. Zhu, J., Zhao, M., Zhang, S., Chang, X.-W., and Precup, D. (2022). "Revis- iting Heterophily For Graph Neural Networks". In: Advances in Neural Information Processing
Systems 35. pp. 1362-1375.
Luan, S., Zhao, M.. Chang, X.WW., and Precup, D. (2019). "Break the Ceiling: Stronger Multi-scale
Deep Graph Convolutional Networks". In: Advances in Neural Information Processing Systems.
Vol. 32. Curran Associates, Inc.
Maurya, S. K., Liu, X., and Murata, (2021). Improving Graph Neural Networks with Simple Architecture Design. arXiv: 2105.07634 [ca, stat].
McCallum, A. K. Nigam, K., Rennie, J., and Seymore, K. (2000). " "Automating the Construction of Internet Portals with Machine Learning". In: Information Retrieval 3.2, pp. 127-163.
Monti, F., Frasca, F., Eynard, D., Mannion, D., and Bronstein, M. M. (2019). Fake News Delection Social Media Using Geometric Deep Learning. arXiv: 1902.06673.
Namata, G. M., London, B., Getoor, L., and Huang, B. (2012). "Query-driven Active Surveying for Collective Classification" In: Workshop on Mining and Learning with Graphs.
Oono, K. and Suzuki, T. (2019). "Graph Neural Networks Exponentially Lose Expressive Power for Node Classification". In: International Conference on Learning Representations.
Paszke, A. et al. (2019). "PyTorch: An Imperative Style, High-Performant
Deep Learning Library". In: Advances in Neural Information Processing Systems. Vol. 32. Curran Associates, Inc.
Pei, H., Wei, B., Chang, K. C.-C., Lei, Y., and Yang, B. (2019). "Geom-GCN: Geometric Graph Convolutional Networks". In: International Conference on Learning Representations. Perko, L. (2001). Differential Equations and Dynamical Systems. Ed. by J. E. Marsden, L. Sirovich, and M. Golubitsky, Vol. 7. Texts in Applied Mathematics. New York, NY: Springer New York.
12

---

Platonov, O., Kuznedelev, D., Diskin, M., Babenko, A., and Prokhorenkova, L. (2023). "A critical look at the evaluation of GNNs under heterophily: Are we really making progress?" In: The Eleventh International Conference
Learning Representations.
Poli, M., Massaroli, S., Park, J., Yamashita, A., Asama, H., and Park, J. (2021). Graph Neural Ordinary Differential Equations. arXiv: 1911. 09522 [cs, stat].
Pozrikidis, C. (2018). The Fractional Laplacian. First. Boca Raton: Taylor & Francis, 2016. I "A
CRC title.": Chapman and Hall/CRC.
Rong, Y., Huang, W. Xu, T., and Huang, J. (2020). "DropEdge: Towards Deep Graph Convolutional
Networks on Node Classification". In: International Conference Learning Representations.
Rozemberczki, B., Allen, C., and Sarkar, R. (2021). "Multi-Scabl Attributed Node Embedding". In: Journal of Complex Networks 9.2.
Rusch, T. K., Chamberlain, B., Rowbottom, J., Mishra, S., and Bronstein, M. (2022). "Graph-Coupled Oscillator Networks". In: Proceedings of the 39th International Conference on Machine Learning.
PMLR. pp. 18888-18909.
Scarselli, F., Gori, M., Tsoi, A. C., Hagenbuchner, M., and Monfardini, G. (2009). "The Graph Neural
Network Model". In: IEEE Transactions on Neural Networks 20.1, pp. 61-80.
Sen, P., Namata, G., Bilgic, M., Getoor, L.., Galligher, B., and Rliassi-Rad, T. (2008). "Collective Classification in Network Data". In: Al Magazine 29.3, 93.
Shi, Y., Huang, Z., Feng. S., Zhong, H., Wang, W., and Sun, Y. (2021). "Masked Label Prediction: Unified Message Passing Model for semi-Supervised Classification". In: Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence. Montreal, Canada: International
Joint Conferences on Artificial Intelligence Organization, pp. 1 1548-1554.
Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R. (2014). "Dropout: Simple Way to Prevent Neural Networks from Overfitting,. In.
Tang, J., Sun, J., Wang. C., and Yang. Z. (2009). "Social Influence Analysis in Large-Scale Networks".
In: Proceedings of the 15th ACM SIGKDD International Conference Knowledge Discovery and
Data Mining. Paris France: ACM, pp. 807-816.
Tong. Z. Liang, Y., Sun, C., Li, X., Rosenblum, D., and Lim, A. (2020). "Digraph Inception
Convolutional Networks". " In: Advances in Neural Information Processing Systems. Vol. 33. Curran
Associates, Inc., pp. 17907-17918.
Tong. Z., Liang, Y. Sun, C., Rosenblum, D. S., and Lim, A. (2020). Directed Graph Convolutional Network. arxi 2004. 13970 [cs, stat].
Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y. (2018). "Graph Attention Networks". In: International Conference on Learning Representations.
Wang. J., Huang, P., Zhao, H., Zhang, Z., Zhao, B., and Lee, D. L. (2018). "Billion-Scale Commodity Embedding for E-commerce Recommendation in Alibaba". In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, KDD '18. New York, NY, USA: Association for Computing Machinery, pp. 839-848.
Wang, X. and Zhang, M. (2022). "How Powerful Are Spectral Graph Neural Networks". In: Proceed- ings of the 39th International Conference on Machine Learning. PMLR. pp. 23341-23362. Wang. Y., Yi, K., Liu, X., Wang, Y. G., and Jin, S. (2022). "ACMP: Allen-Cahn Message Passing with Attractive and Repulsive Forces for Graph Neural Networks". In: The Eleventh International Conference on Learning Representations.
Wu, F., Souza, A., Zhang, T., Fifty, C., Yu, T., and Weinberger, K. (2019). "Simplifying Graph Convolutional Networks". In: Proceedings of the 36th International Conference on Machine Learning. PMLR. pp. 6861-6871.
Wu, X., Chen, Z., Wang, W. W. and Jadbabaie, A. (2022). "A Non-Asymptotic Analysis of Over- smoothing in Graph Neural Networks". In: The Eleventh International Conference on Learning Representations.
Xhonneux, L.-P., Qu, M., and Tang. J. (2020). "Continuous Graph Neural Networks". In: Proceedings of the 37th International Conference on Machine Learning. PMLR, pp. 10432-10441. Xu, K., Hu, W., Leskovec, J., and Jegelka, S. (2018). "How Powerful Are Graph Neural Networks?" In: International Conference on Learning Representations.
Yan, Y., Hashemi, M., Swersky. K., Yang, Y., and Koutra, D. (2022). "Two Sides of the Same Coin: Heterophily and Oversmoothing in Graph Convolutional Neural Networks". In: 2022 IEEE International Conference on Data Mining (ICDM). Orlando, FL, USA: IEEE, pp. 1287-1292. Zhang, X., He, Y., Brugnone, N., Perlmutter, M., and Hirn, M. (2021). "MagNet: A Neural Network for Directed Graphs" In: Advances in Neural Information Processing Systems. Vol. 34. Curran Associates, Inc., pp. 27003-27015
13

---

Zhao, L. and Akoglu, L. (2019). "PairNorm: Tackling Oversmoothing in GNNs". In: International Conference on Learning Representations.
7hu, J., Rossi, R. A., Rao, A., Mai, T., Lipka, N., Ahmed, N. K., and Koutra, D. (2021). "Graph Neural
Networks with Heterophily". In: Proceedings of the AAAI Conference on Artificial Intelligence
32.12, pp. 11168-11176.
Zhu, J., Yan, Y. Zhao, L., Heimann, M., Akoglu, L., and Koutra, D. (2020). "Beyond Homophily
in Graph Neural Networks: Current Limitations and Effective Designs". In: Advances in Neural Information Processing Systems. Vol. 33. Curran Associates, Inc., pp. 7793-7804.
Zou, C., Han, A., Lin, L.,, and Gao, J. (2022). A Simple Yet Effective SVD-GCN for Directed Graphs. arXiv: 2205 09335 [cs].
14

---

Acronyms
AUROC
Area under the ROC curve
DSBM
Directed Stochastic Block Model
FD
Frequency Dominant
FGL
Fractional Graph Laplacian
GAT
Graph Attention Network
GCN
Graph Convolutional Network
GNN
Graph Neural Network
HFD
Highest Frequency Dominant
LCC
Largest Connected Components Lowest Frequency Dominant
LFD
MLP
Multi-Layer Perceptron
ODE
Ordinary Differential Equation
SNA
Symmetrically Normalized Adjacency Singular Value Decomposition
SVD
Notation
Imaginary unit
R(z)
Real part of C
3(2)
Imaginary part of 2 E C
diag(x) (x)
Diagonal matrix with x on the diagonal.
I
Constant vector of all 1s. Transpose of M
MT
M
Conjugate of M
M
Conjugate transpose of M
[M]
Spectral norm of M Frobenius norm of M
M
2
A (M)
Spectrum of M
(M)
Singular values of M
& (x)
Dirichlet energy computed on x
H (g)
Homophily coefficient of the graph G
A B
Kronecker product between A and B
vec (M)
Vector obtained stacking columns of M.
15

---

A Implementation Details
In this section, we give the details on the numerical results Section 6. We begin by describing the exact model.
Model architecture. Let & be directed graphs and X E RNN the node features. Our architec- ture first embeds the input node features via multi-layer perceptron (MLP). We then evolve the features Xg according to a slightly modified version of (3), i.e, x' = -i (t)(t) W for some time t E 00 Ti. In our experiments, we approximate the solution with Explicit Euler scheme with step size h > This leads to the following update rule
X X, ihL W.
The channel mixing matrix is diagonal learnable matrix W C and E R, h C are also learnable parameters. The features at the last time step XI are then fed into second MLP, whose output is used as the final output. Both MLPs use LeakyReLU as non-linearity and dropout (Srivastava et al., 2014). On the contrary, the graph layers do not use any dropout nor non- con-linearity. A sketch of the algorithm is reported in fLode.
Algorithm 1: FLode
% A, Xo are given.
% Preprocessing
Din diag (A1)
2 Doul diag (AT1)
D in A>DD"
U. VH = sVd
% The core of the algorithm is very simple
5 def training_ xxxiinggste((x)::::
Xo input_ttiia. Siput_MLP(xxx (x,,
7
for t E [l....................... T) do
X2 x
UE V
Xy = output Output_MLPPCCC (x7)
10
return XT
Complexity. The computation of the SVD is O(N3). However, one can compute only the first p < N singular values: this cuts down the cost to 2(N2p) p). The memory required store the singular vectors is O(N2), since they are not sparse in general. Each training step has cost of O(N2 K).
Experimental details. Our model is implemented in PyTorch (Paszke et al., 2019), using PyTorch geometric (Fey al., 2019). The computation of the SVD for the fractional Laplacian imple- mented using the library linalg provided by PyTorch. In the case of truncated SVD, we use the function randomized, svd provided by the library extmath from sklearn. The code and instruc- tions to reproduce the experiments are available search. All experiments were run on an internal cluster with NVIDIA GeForce RTX 2080 Ti and NVIDIA TITAN RTX GPUs with 16 and 24 GB of memory, respectively.
GitHub. Hyperparameters were tuned using grid
Training details. All models were trained for 1000 epochs using Adam (Kingma et al., 2015) as optimizer with fixed learning rate. We perform early stopping if the validation metric does not increase for 200 epochs.
A.1 Real-World Graphs
Undirected graphs We conducted 10 repetitions using data splits obtained from (Pei et al., 2019). For each split, 48% the nodes are used for training, 32% for validation and 20% for testing. In all datasets, we considered the largest connected component (LCC). Chameleon, Squirrel, and Film are
16

---

directed graphs; hence, we converted them to undirected. Cora, Citeseer, and Pubmed are already undirected graphs: to these, we added self-loops. We normalized the input node features for all graphs.
As baseline models, we considered the same models as in (Di Giovanni et al., 2023). The results were provided by Pei et al. (2019) and include standard GNNs, such as GAT (Velickovic et al., 2018), GCN (Kipf et al., 2017), and GraphSAGE (Hamilton et al., 2017). We also included models designed to address oversmoothing and heterophilic graphs, such as PairNorm (L. Zhao et al., 2019), GGCN (Yan al., 2022). Geom-GCN (Pei al., 2019), HgGGCN (Zhu, Yan, et al., 2020), GPRGNN (Chien et al., 2021), and Sheaf (Bodnar et al., 2022). Furthermore, we included the graph neural ODE-based approaches. CGNN (Xhonneux al.. 2020) and GRAND (Chamberlain et al.. 2021), as in (Di Giovanni et al., 2023), and the model GRAFF from (Di Giovanni et al., 2023) itself. Finally, we included GREAD (Choi et al., 2023), GraphCON (Rusch et al., 2022), ACMP (Y. Wang et al., 2022) and GCN and GAT equipped with DropEdge (Rong et al., 2020).
Heteriphill-speecfficc Models For ceteroopillyspeccine datasets, we use the same models and results as in (Platonov et al., 2023). As baseline models we considered the topology-agnostic ResNet (He et al., 2016) and two graph-aware modifications: ResNet+SGC(F. Wu et al., 2019) where the initial node features are multiplied by powers of the SNA, and ResNer+adj, where rows of the adjacency matrix are used as additional node features; GCN (Kipf et al., 2017), GraphS (Hamilton et al., 2017);; GAT (Velickovic et al., 2018) and GT (Shi et al., 2021) as well as their modification GAT-sep and GT-sep which separate ego- and neighbor embeddings; H2 GCN (Zhu, Yan, et al., 2020), CPGNN (Zhu, Rossi, et al., 2021), GPRGNN (Chien et al., 202 1). FSGNN (Maurya et al., 2021), GloGNN (X. Li et al., 2022), FAGCN (Bo et al., 2021), GBK-GNN (Du et al., 2022), and JacobiConv (X. Wang et al., 2022).
The exact hyperparameters for FLODE are provided in Table 5.
A.2 Synthetic Directed Graphs
The dataset and code are taken from (Zhang et al., 2021). As baseline models, we considered the ones in (Zhang et al., 2021) for which we report the corresponding results. The baseline models include standard GNNs, such as ChebNet (Defferrard et al., 2016), GCN (Kipf et al., 2017), GraphSAGE (Hamilton et al., 2017), APPNP (Gasteiger et al., 2018), GIN (Xu et al., 2018 8), GAT (Velickovic et al., 2018), but also models specifically designed for directed graphs, such as DGCN (Tong, Liang, Sun, Rosenblum, et al., 2020), DiGraph and DiGraphIB (Tong, Liang, Sun, X. Li, et al., 2020), MagNet (Zhang et al., 2021)).
The DSBM dataset. The directed stochastic block model (DSBM) is described in detail in (Zhang et al., 2021, Section 5.1.1). To be self-contained, we include short explanation.
The DSBM model is defined as follows. There are N vertices, which are divided into ne clusters (C1, Cp.),..,, each having an equal number of vertices. An interaction is defined between any two distinct vertices, u and v, based on two sets of probabilities: {a and {
The set of probabilities {as.y) is used to create an undirected edge between any two vertices u and v, where u belongs to cluster C; and v belongs to cluster Cj. The key property of this probability set is that os = aj.i. which means the chance of forming an edge between two clusters is the same in either direction.
The set of probabilities {3 } is used to assign direction to the undirected edges. For all i,j
{I...., J. we assume that B + 3, = 1 holds. Then, to the undirected edge (u, v) is assigned the direction from and the direction from v to u with probability 83..-
with probability Br,j if belongs cluster C, and belongs to cluster Cy.
The primary objective here is to classify the vertices based on their respective clusters.
There are several scenarios designed to test different aspects of the baseline models and our model. In the experiments, the total number of nodes is fixed at N = 2500 and the number of clusters is fixed at n. = 5. In experiments, the training set contains 20 nodes per cluster, 500 nodes for validation, and the rest for testing. The results are averaged over 5 different seeds and splits.
17

---

Table Test accuracy on node classification: top three models indicated as
SE
2nd 3rd
(a) Undirected graphs.
Squirrel
Chameleon
Citeseer
Cora
GGON
37.54
1.56
55.17 + 1.58
71.14
1.84
77.14 + 1.45
80.
0.37
87.95 1.05
GPRGNN
31.61 1.24
46.58
1.71
77.13
1.67
87.54
57.95
FAGCN
2.15
77.11
1.57
37.44
1.30
38.47 + 1.58
63.86
3.04
77.33 + 1.48
90.15 + 0.43
[S8.37 1.26
Goom-GCN
31.59 1.15
38.15 1: 0.92
2.81 60.00
78.02 1 1.15
89.95
0.47
85.35 1 1.57
PairNorm
27.40 1.24
50.44 1: 2.04
62.74
73.59 1: 1.47
87.53 1: 0.44
85.79 1: 1.01
GraphSAGE
34.23 = 0.99
41.61 0.74
18.73
76.01 + 1.30
88. + 0.50
86.90 = 1.04
GCN
27.32 = 1.10
53.43 2.01
64.82 2:24
76.50 + 1.36
88.42 4 0.50
86.98 = 1.27
GAT
27.44 = 0.89
40.72 1.55
60.26 2.50
76.55 + 1.23
87.30 + 1.10
86.33 = 0.48
MLP
36.03 0.70
28.77 = 1.56
46.21 2.99
74.02 - 1.90
75.69 + 2.00
87.16 = 0.37
CGNN
35.95 0.86
29.24 1: 1.09
46.89
1.66
76.91 1. 1.81
87.70 0.49
87.10 1.35
GRAND
35.62 1.01
40.05 1: 1.50
54.67 1. 2.54
76.46 1. 1.77
89.02 0.51
87.36 0.96
Sheaf (max)
37.81
56.31 1.32
76.70 1.57
GRAFFNL
35.96
59.01 + 1.31
1.47 71.38
76.81 +: 1.12
0.50
87.81
GREAD
37.90 1.17
[59.22 1.44
11.38
1.30
77.60 1.81
00.23 0.55
88.57
0.66
GraphOON
35.58 1.24
35.51 + 1.40
49.63
1.89
76.36
2.67
88.01 0.47
87.22
1.48
ACMP
34.93
1: 1.53
57.39
2.09
76.71
87.79
87.71
GCN+DropEide
1.77
2.04 59.06
76.37
83.54
Gro+Dre
28.95
0.76
41.27 + 1.76
58.95
2.13
76.13 + 2.20
86.91 0.45
83.54 1.06
FLODE
37.16 1.42
6:123 1.8
73,60 1.55
1.62
86.11 1.17
(b) Directed graphs.
Film
Squirrel
Chameleon
ACM
36.89 +: 1.18
54.4 +: 1.88
67.08 :: 2.04
HLP
34.59 + 1.32
74.17 1.83
77 77.48 = 1.50
ESGNN
35.67 1 0.69
73.48 2.13
78.14 1.25
GRAFF
37.11 1: 1.08
58.72
0.84
71.08 :: 1.75
FLODE
37.41 + 1.06
74.03 + 1.58]
77.98 + 1.05
(c) Heterophily- graphs Minesweeper, Tolokers and Questions the evaluation metric is AUROC.
Roman-empire
Minesweeper
Tolokers
Questions
ResNet
63.88 +: 0.38
1.39 70.88 0.90
50.89
72.95 +: 1.06
70.34 0.76
ResNeH4SGC
73.90 1: 0.51
80.70 1 0.97
75.81 0.96
RosNot+adj
52.25 + 0.10
50.12 + 0.83
78.78 + 1.11
75.77 = 1.21
73.69 + 0.74
89.75 0.52
83.64 + 0.67
76.09 1.27
GraphSAGE
85.74 + 0.67
93.51 0.57
82.43 - 0.44
76.44 U.62
GAT
80.87
+ 0.30
92.01 + 0.68
83.70 + 0.47
77.43 + 1.20
GAT-sep
88.75 + 0.41
93.91 1: 0.35
83.78 1 0.43
76.79 : U.71
GT
86.51 + 0.73
91.85 1: 0.76
83.23 1 0.64
77.95 0.68
GT-sep
87.32 + 0.39
0.47 92.29
82.52 = 0.92
78.05 0.93
FAGCN
60.11 + 0.52
89.71 0.31
73.35 - 1.01
63.59 = 1.46
CPGNN
63.96 + 0.62
52.03 5.46
73.36 - 1.01
65.96 = 1.95
HgGCN
64.85 + 0.27
86.24 2 0.61
72.94 + 0.97
55.48 0.91
ESGNN
79.92 + 0.56
90.08 01.88:00 0.70
82.76 + 0.61
78.86 0.92
GkoGNN
59.63 + 0.69
51.08 1: 1.23
73.39 1 1.17
65.74 1.19
FAGON
65.22 + 0.56
0.73
77.75 1.05
77.24
1.26
GBK-GNN
74.57 0.47
81.01 0.67
74.47
IncobiConv
0.42
89.66
0.65
73.88
FLODE
74.97 0.53
92.43 0.51
84.17 1 0.58
78.39 1.22
18

---

Table 5: Selected hyperparameters, learned exponent, step size, and Dirichlet energy in the last layer for real-world datasets.
(a) Undirected.
Dataset
Film
Squirrel
Chameleon
Citeseer
Pubmed
Cora
learning rale
10-3
2.5 10-3
5 1-100
10-2
10-2
10-2
weight decay
5 10
5. 10
10
10
10
5 10
channels
layers
encoder layers
1
docoder layors
2
input dropout
1.5 10
0
5 10-2
decoder dropout
10
10 1
0
10 11
exponent
0.003
0.17 = 0.03
0.15
0.92 0.03
0.82
step size
0.991
0.002
1.08 0.01
1.22 0.03
1.04 0.02
1.12 0.02
1.06 0.01
Dirichlet energy
0.246
0.006
0.40 0.02
0.13 0.03
0.021 1: 0.001
0.015 0.001
0.0227 0.0006
(b) Directed.
Dataset
Film
Squirrel
Chameleon
leamingg rate
10
2.5 10
10
weight decay
5 10-1
10-1
10-3
hidden channels
64
64
layers
6
layers layers
input dropout
10
decoder dropout
0.1
10-1
exponent
1.001
0.28 + 0.06
step size
0.990 0.002
1.22 = 0.02
1.22 = 0.05
Dirichlet energy
0.316 0.005
0.38 1 0.02
0.27 1 0.04
(c) Heterophily-speciffic graphs
Roman-empire
Minesweeper
Tolokers
Questions
learning rate
10-3
10-3
10-3
10-2
weight decay
0
0
5 10
hidden chamnels
512
512
512
128
layers
4
5
eucoder layers
decoder layers
input dropout
decoder dropout
0
0
expoment
1: 0.046
step
1:
Dirichlet energy
0,059 0.003
0.173 = 0.019
0.135 4 0.013
0.092 + 0.039
19

---

Table 6: Node classification accuracy of ordered DSBM graphs: top three models as (a) Varying edge density.
2 and 3rd
0.1
0.08
0.05
ChebNet
19.9 + 0.6
20.0 + 0.7
20.0 0.7 0.7
GCN-D
68.9 = 2.1
67.6 + 2.7
58.5
APPNP-D
97.7 1.7
95.9 + 2.2
90.3
GraphSAGE-D
20.1 1.1 19.9
0.8
19.9
GIN-D
57.3 1: 5.8
55.4
1 5.5
50.9 7.7
GAT-D
42.1 + 5.3
39.0
+ 7.0
37.2
- 7.2
+
64.4
DiGrrphh
82.1
77.7 + 1.6
66.1
2.
DiGmahIB
99.2 11200.1 0.5
97.7 0.7
89.3
1.7
MagNet
99.6 0000002 0.2
98.3
0.8
91. 1.2
FLODE
99.3 = 0.1
98.8
(b) Varying net flow.
0.05
0.10
0.15
0.20
0.25
0.30
0.35
0.40
ChebNet
19.9 0.7
20.1 +
0.6
20.0 0.6
20.1 + 0.8
19.9
0.9
20.0 +
0.5
19.7 +
0.9
20.0 + 0.5
GCN-D
68.6 2.2
74.1
7.
74.9 1.3
72.0
0..41.0 1.6
58.1
45.6
APPNP-D
97.1 1.8
2.4
7.88 3.6
69.4
3.9
59.6
4.9
51. 4.5
39.1 5.3
GraphSAGE-D
20.2 1.2
2..0
20.0 2.80
19.6
19.8 0.7
19.9 0.9
19.9 0.8
GIN-D
57.9 1 6.3 48.0
11.4
32.7
+ 12.9
26.5
10.0
23.8
2
6.0
20.6
3.0
20.5
2.8
19.8
0.5
GAT-D
42.0
32.7
5.1
25.6 3.8
19.9 = 1.4
20.0
1.0
19.8
+ 0.8
19.6
0.2
19.5
0.2
DGCN
81.4 1.1
8.77
85.5 1.0
86.2 = 0.8
84.2 111
78.4 13
09.0 1.6
6IKE 160
DiGraph
82.5 1.4
1.9
81.9
79.7
1 1.3
73.5 1: 1.9
67.4 1 2.8
57.8 + 1.6
43.0 7.1
DiGraphIB
99.2
07.9
94.1+
88.7
2.0
82.3 = 2.7
70.0
+ 2.2
57.8 6.4
41.0
MagNet
99.6 0.2
99.0
100010 1.0
97.5
912 1.6
88.7
199
79.1
2.9
68.8
2.1
518
FLODE
99.3 0.1
0.1
98.3 U.1
96.7 + 0.2
92.8 02820.1 0.1
87.2 0.3
77.1 0.5
0.5
63.8 0.3
50.1 0.5
Following Zhang al. (2021), train model in both experiments for 3000 epochs and use carly-stopping if the validation accuracy does not increase for 500 epochs. We select the best model based on the validation accuracy after sweeping over few hyperparameters. We give exact numerical values for the experiments with the standard error in Table 6a and refer to Appendix A.2 for the chosen hyperparameters.
DSBM with varying edge density. In the first experiment, the model is evaluated based on its performance on the DSBM with varying as a', a E {0.1, 0.08, 0.05} for + j. which essentially changes the density of edges between different clusters. The other probabilities are fixed at a; 0.5, Br, = 0.5 and Br, = 0.05 for j. The results are shown Figure 6 with exact numerical values in Table 6a.
DSBM with varying net flow. In the other scenario, the model is tested on how it performs when the net flow from one cluster to another varies. This achieved by keeping a,, = 0. constant for all i and j. and allowing 3rg to vary from 0.05 to 0.4. The other probabilities are fixed at a,,i = 0.5 and Br, 0.5. The results are shown in Figure 6 with exact numerical values in Table 6b.
A.3 Ablation Study
We perform an ablation study on Chameleon and Squirrel (directed, heterophilic), and Citeseer (undi- rected, homophilic). For this, we sweep over different model options using the same hyperparameters
20

---

Table Selected hyperparameters for DSBM dataset.
(a) Varying edge density.
0.1
0.08
0.05
leurning rale
5 10
5 5-10-3
5 10-3
decay
1- 10
1-10
5 10
input dropout
1 1-10-1
2-10-1
1 1.10-1
decoder dipour
10
10
hidden channels
256
256
256
(b) Varying net flow.
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
learning rate
5
10-3
10
1
10
10
docay
1 1-10-3
1-
1
1.
5
input dropout
1 -00
1- 10
2 2-10-
1 10
1. 10
2 10
5 10
10
decoder dropout
1 1.10
1. 10
10
5 10
5 10
10
2 10
hidden channels
256
256
256
256
256
256
256
256
via grid search. The test accuracy corresponding to the hyperparameters that yielded maximum validation accuracy is reported in Table 8.
The ablation study Chameleon demonstrates that all the components of the model (learnable exponent, ODE framework with the Schrödinger equation, and directionality via the SNA) contribute to the performance of FLODE. The fact that performance drops when any of these components not used suggests that they all play crucial roles in the model's ability to capture the structure and evolution of heterophilic graphs. It is important to note that the performance appears be more dependent on the adjustable fraction in the FGL than on the use of the ODE framework, illustrating that the fractional Laplacian alone can effectively capture long-range dependencies. However, when the ODE framework is additionally employed, noticeable decrease in variance is observed.
From Theory Practice. We conduct an ablation study investigate the role of depth on Chameleon, Citeseer, Cora, and Squirrel datasets. The results, depicted in Figure 8, demonstrate that the neural ODE framework enables GNNs scale to large depths (256 layers). Moreover, we see that the fractional Laplacian improves over the standard Laplacian in the heterophilic graphs which is supported by our claims in Section 5.2. We highlight that using only the fractional Laplacian without the neural ODE framework oftentimes outperforms the standard Laplacian with the neural ODE framework. This indicates the importance of the long-range connections built by the fractional Laplacian.
We further demonstrate the close alignment of our theoretical and experimental results, which enables us to precisely anticipate when the models will exhibit HFD or LFD behaviors. In this context, caleulate parameters (according Theorem D.5) and illustrate each depth the expected and observed behaviors. For Squirrel and Chameleon, which are heterophilic graphs, we observe that both their theoretical and empirical behaviors are HFD. Additionally, the learned exponent is small. In contrast, for Cora and Citeseer, we see the opposite.
Finally, employ the best hyperparameters in Table 5a solve both fractional heat and Schrödinger graph ODEs, further substantiating the intimate link between our theoretical advancements and practical applications.
21

---

Table Ablation study node classification task: top models are indicated (a) Chameleon (directed, heterophilic).
1st and 2nd
Update Rule
Test Accuracy
Dirichlet Energy
X(+1 X ihL"x, X; W
77.79 + 1.42
0.213 (t=5)
D X X ihL X; W
75.72 + 1.13
0.169 (t=6)
X(+1 -i L W
77.35 2.22
0.177 (t=4)
Xx+1 = -i X W
69.61 1.59
0.178 (t=4)
Xt+1 X i X, W
73.60 1.68
0.131 (t=4)
X(+1 = X ihL X W
70.15 0.86
0.035 (t=4)
X,+1 = i L'x, W
71.25 - 3.04
0.118 (t=4)
X(+1 i L X W
67.19 + 2.49
0.040 (t=4)
Xx+1 Xx hL W
77.33 + 1.47
0.378 (t=6)
Xx+1 X hL X, W
73.55 0.94
0.165 (1=6)
D
X(+I L W
74.12 3.60
0.182 (t=4)
Xe+1 = L X W
68.47 + 2.77
0.208 (1=4)
(b) Squirrel (directed, heterophilic).
Update Rule
Test Accuracy
Dirichlet Energy
X(+1 = X; ihL'x, W
74.03 1.58
0.38 0.02
D X(+1 = Xx ihL X W
64.04 + 2.25
0.35 + 0.02
X(+1 i L W
64.25 + 1.85
0.46 0.01
Xx+1 i L X W
42.04 1.58
0.29 0.05
Xx+1 X, ihL*x x, W
64.23 1.84
0.40 0.02
X(+1 X ihL X, W
55.19 1.52
0.26 0.03
X,+1 = i L'x, W
61.40 2.15
0.43 0.01
X(41 = i L X, W
41.19 + 1.95
0.20 - 0.02
X(41 = Xx hL"x, W
71.86 + 1.65
0.50 + 0.01
X(+1 Xt hL x, W
59.34 + 1.78
0.43 - 0.03
D
X(+1
L W
42.91 7.86
0.32 0.08
Xx+1 =
X,
35.37 1.69
0.25 0.05
X(+1 X, IL x, W
62.95 2.02
0.61 0.08
X(+1 = X2 h X; W
52.19 1.17
0.51 0.07
X(+1 = L W
59.04 + 0.02
0.44 - 0.02
X(+1 L X; W
39.69 + 1.54
0.20 0.02
(c) Citeseer (undirected, homphilic).
Update Rule
Test Accuracy
Dirichlet Energy
Xc+1 X₁ ihL x, W
78.07 + 1.62
0.021 (t=5)
X(+1 X, ihL x, W
77.97 + 2.29
0.019 (t=4)
X(+1 = -i L
77.27 2.10
0.011 (t=6)
Xt+1= = -i L X, W
77.97 +2.23
0.019 (t=4)
22

---

Table 9: Learned and spectrum of W. According Theorem 5.3. denote FD == xx (W) (A1 (W) and FD = 3(A((W)) (Ar
(W)) for the frac- tional heat (H) and Schrödinger (S) graph ODEs, respectively. The heterophilic graphs Squirrel and Chameleon exhibit HFD since FD O, while the homophilic Cora, Citeseer, Pubmed exhibit LFD
since FD O.
Squierel -0.8896
Chameleon -0.9337
(L
0.9488
-0.5022
-0.6537
-0.4826
1.008 + 0.007
0.19 0.0G
0.37 + 0.14
0.89 + 0.06
1.10 0.08
0.89 + 0.01
(r(W)
-2.774
0.004
-1.62 0.03
-1.81 4: 0.02
-1.76 0.01
-1.68 0.08
-1.81 0.01
H
(W)
2.858 0.009
2.32
FD
0.367 0.001
-0.54 0.02
-0.42 2: 0.04
0.52 + 0.02
0.97 0.09
0.60 0.01
((A)
-2.95 0.001
-1.70 +: 0.04
1.74
-1.78 0.01
(W))
2,880 1: 0.002
2.20
2.46 1: 0.02
2.29 +: 0.07
0.98 0.09
2.30 1: 0.02
0.4945 0.0001
-0.48
-0.62 0.03
0.46 = 0.06
1.03 0.05
0.59 0.01
X1+1 X
/L
W
X =
X
WLX:
X = L'
X141 Lx W
Cora
Test Test accuracy
0.8
Citeseer
0.6
Chameleon
0.4
Squirrel
0.2
10
X0
100
/8
10
1-00
LED LFF
LED HED
0.5
Figure 8: Ablation study the effect of different update rules and different number of layers undirected datasets. The x-axis shows the number layers 2 for L E (O....,8) 88. FD calculated according to Theorem 5.3.
23

---

B Appendix for Section
Proposition 3.3. Let 9 be directed graph with SNA L. For every X E X(L), holds [|] S and = -L)(
Proof. We show that the numerical range W(L) Lx*Lx x*x = 1) satisfies W(L) C [[-1, 1]. As (L) contains all eigenvalues of L the thesis follows.
Let A be the adjacency matrix of G and x E CV with x"x= 1. Applying the Cauchy-Schie Souchy-Schhwrt inequality in (2) and (3), we get
N
N
x
|zy|
x Lx
Lx
d: doul
t=1
j=1
22
as.
ague
==1
N
|zz 1 dout
Oi,
aij
j=1
رzح dour
aid
j=1
(3)
zzy 2 dgu
asd
Σ
where we used oř.5 = aigz We have -i=1 |x:|² 2- xxx such that W(L) [-1, 1] follows. The Av A)v.
second claim follows directly by (I- (L-L)
Proposition 3.5. Let G be a directed graph with SNA Then E X(L) if and only if the graph is weakly balanced. Suppose the graph is strongly connected; then - E X(L) if and only if the graph is weakly balanced with an even period.
Proof. Since the numerical range is only a superset of the set of eigenvalues, we cannot simply consider when the inequalities (1) (3) in the previous proof are actual equalities. Therefore, have to find another way to prove the statement. Suppose that the graph is weakly balanced, then
N
ky
hy d
Σ
aij
O, E (1,..,
jout
j=1
We will prove that k = (
is an eigenvector corresponding to the eigenvalue 1,
(Lk);
Gi,j
Ki
k
1
Σ ais
k; ki
din
For the other direction, suppose that there exists X E R such that * 0 and x i E {1....,N) NN
Lx. Then for all
N
0 (Lx); 24
Cijj
24
Di,j Fy
Dijj in
din
d
d d
j=1
24

---

N
hence, the graph is weakly balanced.
By Ferron-Frobenius theorem for irreducible non-negative matrices, one gets that L has exactly h eigenvalues with maximal modulus corresponding to the roots of the unity, where the period of L. Hence, -1 is eigenvalue of L if and only if the graph is weakly balanced and h is even. G
Proposition 3.6. For every x E C XX, we have
R (trace (x# (I - L))
aj
your j
2,j=1
Moreover, there exists + 0 such that ((x) 0 if and only if the graph is weakly balanced.
Proof. By direct computation, holds
N
I
Oi
2
d
pour
K
Ti,k
J...
Di.j
2
di
1,
N
K
Σ
Ti,k
j.jk
J.
O5.j
2
ddd
your
1-=1
|20.A 2
K
O5.j
+
diTh
1=
1,==1 k=1
I' j.
Σ
Σ
04,
d dou
di dzue
W=
N
K
1
N
K
T
N
K
It 2
Fi.h
2
+
2
Gijj
aij
2
2
2
du dgua
/ di dgu
i=1
j=l
k=1
k=
b R=
N
K
N
K
N
K
J. (x. )kss
Ot.j
2
d dgu
V /d ar
WJ=1
k=
1==1k=1 k=1
K
N
(xtt)
N
K
(x+t) (. / din drm
as.j
d dou
1,==1
1.f=1
K
N
K
R
2
Oi,j
ddf dou
k=
1,==
km
Rt (trace (xH (I L)))
The last claim can be proved as follows. For simplicity, suppose x E RNN, The " = is clear since one can choose x to be k. To prove the = " we reason by contradiction. Suppose there exists x f 0 such that (x) = 0 and the underlying graph is not weakly connected, i.c..
N
T
VX # 0,
aijj
> 0, Vi [)....,N,,
mm
i=1
25

---

Then, since X + 0,
1 N
0=
6x)
Σ
32
2j
4
your
1.
1
Σ
24
22
Σ ord
ai
d
Vd
V d
N
NV
Σ 1=1
36g
V d
1=1
N
Σ
Jj
aij
Vd
V
0,
where we used Cauchy- Schwartz and triangle inequalities.
We give the following simple corollary. Corollary B.1. For every x RRN it holds 6(x) = 1 (vec(x)) (I L))vec(x)))
C Appendix for Section 4
In this section, provide some properties about FGLs. The first statement shows that the FGL of normal SNA only changes the magnitude of the eigenvalues of L.
Lemma C.1. Ler M be a normal matrix with eigenvalues .I,..., AN and corresponding eigenvectors
Vis..ry VN. Suppose M LER" is its singular value decomposition. Then it holds
Σ = Al L V, R = Vexp (1e) θ diag ((0)). 0g = atan2 (RA,, 31x)
Proof. By hypothesis, there exist a unitary matrix V such that M VA VH then
M VA VA' V A VH
M REL L
Therefore, Σ A and V
M R V
Finally, we note that it must hold R = V exxp(ii) where θ diag ((atan2(RA,, (Aa)) and atan2 is the 2-argument arctangent.
We proceed by proving Theorem 4.1, which follows the proof of a similar result given in (Benzi al., 2020) for the fractional Laplacian defined in the spectral domain of an in-degree normalized graph Laplacian. However, our result also holds for directed graphs and in particular for fractional Laplacians that are defined via the SVD of a graph SNA.
Lemma C.2. Let M E Raxa with singular values (MM) C [a, b). For f : (a, b] - RR define f(M VH where M UEVH is the singular value decomposition of M. If f has
modulus of continuity us and d(i,j) 2 2, it holds
7² [f(M) - 1+
(d(i,)) -(0
w
2
2
26

---

Proof. Let [a, b] - FR be any function, then
(l(() g(M)l2 ||/ Ug Ug(E)V"Il
= -If(2²) -(() =11((2) =())lloc,
The second equation holds since the 2-norm is invariant under unitary transformations. By Jackson's s Theorem, there exists for every m 2 polynomial Pm of order m such that
(I(M) Pm (M) S || Pmmlloo (a. S 1
b- 2
2
2m
ij E (l,...,m.) If J(i,)) m + 1, then any power M up to order m has a zero entry in (i,f), i.e., (M").. 0. Hence, / (M) [.j - f(M). Pm (M)i. and we get
(f S [lf(M) - 2(M)ll2 S w 1+ n² _
b (b-a)
π
a d( d( 1111
= 1
2m
2
2
from which the thesis follows.
Finally, we give proof of Theorem 4.1, which is a consequence of the previous statement.
Proof of Theorem 4.1. The eigenvalues of L are in the unit circle, i.e., |L]| S Hence, [LL*| S 1, and the singular values of L are in [0, 1]. By Lemma C.2 and the fact that f(x) x has modulus of continuity (t) 1a the thesis follows.
D Appendix for Section 5
In this section, we provide the appendix for Section 5. We begin by analyzing the solution of linear matrix ODEs. For this, M CCNN, For Xp E CN consider the initial value problem
(1 = -Mx(t), x(0) = Xo-
(5)
Theorem D.1 (Existence and uniqueness of linear ODE solution). The initial value problem given by (5) has a unique solution x(t) E C for any initial condition Xo C
The solution of (5) can be expressed using matrix exponentials, even if M is not symmetric. The matrix exponential defined as:
(-M)
exp(
Mt)
h!
where M is the k-th power of the matrix M. The solution of (5) can then be written as
x(1) = expp(MM)xo.
(6)
D.1 Appendix for Section 5.1
In this section, we analyze the solution to (2) and (3). We further provide proof for Theorem 5.3. We
begin by considering the solution the fractional heat equation (2). The analysis for the Schrödinger equation (3) follows analogously.
The fractional heat equation x'() Lx W can be vectorized and rewritten via the Kronecker product as
vec(x) W (e)e
(7)
In the undirected case L and I are both symmetric, and the eigenvalues satisfy the relation (I - L) the same for L and (L). The corresponding eigenvectors thx and vA L) be chosen to be L. In the following, we assume that these eigenvectors are orthonormalized.
If L is symmetric, (sh (L), 1/ (L) is an orthogonal matrix containing the eigenvectors of L, and is the diagonal matrix of cigenvalues.
can decompose it via the spectral theorem into L UDUT, where U =
27

---

Due to Lemma C.I, the fractional Laplacian L can be written as L = U a (D)U", where f E the map + sign(z) |x and is applied element-wise. Clearly, the eigendecom- and the corresponding
position of L is given by the eigenvalues (A))) (A)((L)) (aN eigenvectors
Now, by well-known properties of the Kronecker product, one can write the cigendecomposition of W & L as
(), fa
lEft........ (v,(W)
(1 (E(......)
Note that A (L) and, since trace (L) O, the SNA has at least negative eigenvalue. This property useful since it allows to retrieve of the indices (r, 1) corresponding to eigenvalues with minimal real (or imaginary) parts in simple way.
The initial condition vec(xo) can be decomposed as
voc(xo) Σ
Gr v (W) & (L), Cr,1 Gr (vec(Xo), (h) & (h()W)
Then, the solution (tc(x)(t) of (7) can be written as
vec(x)(1) ΣΣ Cr exp (-tòr (W) (Ar (L))) (WW)) & (L))
The following result shows the relationship between the frequencies of L and the Dirichlet energy and serves as a basis for the following proofs.
Lemma D.2. Let G be a graph with SNA L. Consider x(t) E N such that there exists lp NNK \ [0] with
vec(x (1)
1-000 vec(4),
[f)
and @ (--))vec(()) = Avec((p). Then,
x(t)
R((A)
8
2-+00
[x(()]l2
2
Proof. As vec((p) is the limit of unit vectors, vec(yp) is unit vector itself. We calculate its Dirichlet energy.
1
((vec((p) = (( 60 (I L))vec(p))
2 R R(A (evc())) (A)
Since + ((x) is continuous, the thesis follows.
Another useful result that will be extensively used in proving Theorem 5.3 is presented next.
Lemma D.3. Suppose x(t) can be expressed as
x(1)=
Ck. exp ( (tt V& & Wn.
R=1n=1
for some choice of Ck,w d {vx). ). {w,). Let (a, b) be the unique index of dt. with minimal real part and corresponding non-null coefficient Ck i.e.
(a,b) = arg min (kn))[[[x[N]
(R (A CR. # 0}
Then
x(t)
100 Cab 2 W.
xx(t)ll2
Ica. V. & wall2
28

---

Proof.
The key insight is to separate the addend with index (a, 6). It holds
N
x(t) Chem exp (-t Vn @ Wm
= exp(-te -a,b))
Ca.bVa 80 Wb +
Σ
Ckaa exp (-6 (AR,D Aab)) v Wn
(kx))e/K]x[N) (M.n))((a,b)
We note that
lim eep (-t (ARm )a,b)))
lim lexp (-t (Ak, )a,b)) exp (-it (Ak, )a,b))
lim exp (-t (Akm )a,b))
0,
for all (k, # (a, b), since R (AR,
)a,b) O. Therefore, one gets
x(t)
t0000 Ca.b Va @ W;
x(11ll2
|Ca.6 Va @ wall2
where the normalization removes the dependency exp (-1 Aa.b)
When Aa,b not unique, it is still possible derive convergence result. this case, x will converge to an element in the span generated by vectors corresponding to λa i.e.
Σ
Ca.h Va @ Wb
x(t)
1-+00
(a.)EA
[x(1)l1
Σ
Ca.b Va @ woll2
(a.))Ed
where == = ((k, R(AA..)) = R(A CA,n * 0}.
A similar result to Lemma D.3 holds for slightly different representation of x(t). Lemma D.4. Suppose xl) can be expressed as
K N
x(t) =
ΣD
Σ
Ck.. exp (it VA W..
k=
for some choice of Ck.n. A {vR}. {wn). Let (a, b) be the unique index of An with minimal imaginary part and corresponding non-null coefficient Ck i.e.
(a, = arg min (AE,))) Ck # 0} (B.))E[K]([N]
Then
x(t)
10000
Ca.b Va & Wt
[x(+)]la
||ca,b Vo @ woll2
Proof. The proof follows the same reasoning as in the proof of Lemma D.3. The difference that the dominating frequency is the one with the minimal imaginary part, since
R (i Ak.) - (Ak,))
and, consequently,
arg max [K]])[K]x[N]]
(R (i Ak,n))
arg min (bx))E[[K|x())
(S (Ak,a))
29

---

D.1.1 Proof of Theorem 5.3
We denote the eigenvalues of L closest to 0 from above and below as
(L)) arg min (A<(L) AL 0}
(9)
A- (L) = arg max (L
(L) < 0}
We assume that the channel mixing W IR kx and the graph Laplacians L, L E NNN are real matrices. Finally, we suppose the eigenvalues of a generie matrix M are sorted in ascending order, i.e., (M) & AyMM), i < j.
We now reformulate Theorem 5.3 for the fractional heat equation (2) and provide its full proof, which follows similar frequency analysis to the one (Di Giovanni et al., 2023, Theorem B.3)
Theorem D.5. Let G be an undirected graph with SNA L. Consider the initial value problem in (2) with channel mixing matrix W E R and a E R. Then, for almost all initial conditions
XD N the following is satisfied.
(a 0) The solution to (2) HFD if
fa < <>i(W),
and LFD otherwise.
(a < 0) The solution to (2) is (1 (L))--DD if
(x (A- (L)) (A,
and (1 (1-Aa 11- otherwise.
Proof of (a > 0). As derived in (8), the solution of (2) with initial condition Xy can be written in a vectorized form as
vec(x)(t) exp (-t (-t L") vec(xo)
ΣΣ Cra exp (- (W) fo (A((L)) d @ (L().
where X, (W) are the eigenvalues of W with corresponding eigenvectors vl-(W) and (L) are the eigenvalues of with corresponding eigenvectors (L.) The coefficients Cr22 are the Fourier coefficients of xo, i.c.,
Cr == (vec(xo), v. (W) @ (L))
The key insight is to separate the eigenprojection corresponding to the most negative frequency. By Lemma D.3, this frequency component dominates for t going to infinity.
Suppose
(x fo (NT(L) < (W) (Ax((L))))) (L) (W)
In this case, Ax (W) fa (A1 (L)) is the most negative frequency. Assume for simplicity that Ax (W) has multiplicity one: the argument can be applied even if this is not the case, since the corresponding eigenvectors are orthogonal for higher multiplicities.
For almost all initial conditions Xo, the coefficient CK,1 is not null; hence
vec ())
t0000 CK,1 t (W)
@ th (L)
vec (1)
[CK,1 VK (W) @ th (L) 2
By standard properties of the Kronecker product, we have
(I & L) (W (W) @ th (L)) = WI (W)) @ th (L) (L) th (W) 2 v (L)
(10)
ie., dK (W) @ 2/1 (L) is eigenvector of & L corresponding
the eigenvalue X1 (L). Then,
by Proposition 3.3, , (W) @ wi (L) is also an eigenvector of & I L corresponding to the eigenvalue 1 (L) A [ L). An application of Lemma D.2 finishes the proof.
Similarly, we can show that if a > and Ax (W) fa (A1 (L)) X (W) the lowest frequency component A (I- L) is dominant.
30

---

Proof of (a < 0). In this case either fo (x+ (L)) A1 (W) or (A_ (L)) Ag (W) are the most neg
ative frequency components. Hence, (A_ (L)) A (W) > > (A_ (L) X (W) the frequency
la (X+ (L)) A (W) is dominating and otherwise the frequency la (X- (L)) Ax (W). We can see
this by following the exact same reasoning of (i).
Remark D.6. In the proof of (ax < <0), we are lacitly assuming that L has only non-zero eigenvalues. not, we truncate the SVD and remove all zeros singular values (which correspond zeros eigenvalues). In doing so, we obtain the best invertible approximation of L to which the theorem can be applied.
We now generalize the previous result to all directed graphs with normal SNA.
Theorem D.7. Let G be strongly connected directed graph with normal SNA L such that (L) R. Consider the initial value problem in (2) with channel mixing matrix W RRKK and
O. Then, for almost all initial values TR. N the solution to (2) is HFD if
Axc (L))) < (W))
and LFD otherwise.
Proof. Any normal matrix is unitary diagonalizable, i.c., there exist eigenvalues AI...., Av and corresponding eigenvectors Vi...., Vy such that VAVH. Then, by Lemma C.1, the singular value decomposition of L is given by L UEVH, where
Σ = A U
V exp (ie) θ diag (0.)) 0: = alan2 (RAi, 3A.)
Hence,
= -E VIA exp (ie) VH.
Then, equivalent to the derivation of (8), the solution to the vectorized fractional heat equation
W @ ve)e
given by
K Σ
vec vec(x) (x) =
Crd exp (-1Ar (W) (Ar (L))) (W))) 8 (L.)
1=1
with
(A) = |A(L).| exp(i0)).
Now, equivalent to the proof of Theorem 5.3, we apply Lemma D.3. Therefore, the dominating frequency is given by the eigenvalue of W @ L with the most negative real part. The eigenvalues of W @ La are given by X, (W) (L))) for r 1,..., K,I .....,, N. The corresponding real parts given by
RAAr (W) I. (A)))) = A, (W) a cos(O.) = A, (W) -- | R(((L);))
By Perron-Frobenius, the eigenvalue with the largest eigenvalues is given by AL E R. Hence, for all | = 1,...,.1 N,
I cos(0)) < |A(L)n1
Similarly, for all / = 1,..., N with R(((L))) < 0,
| c0: S A(
Thus, the frequency with the most negative real part is either given by (W) f, (A1 (L)) or X f, (L)). The remainder of the proof is analogous to the proof of Theorem D.7.
In the following, we provide the complete statement and proof for the claims made in Theorem 5.3 when the underlying ODE is the Schrödinger equation presented in (3).
31

---

Theorem D.8. Let 9 be undirected graph with SNA L. Consider the initial value problem in (3) with channel mixing matrix W E C and R. Suppose that W has least one eigenvalue with non-zero imaginary part and sort the eigenvalues of W in ascending order with respect to their imaginary part. Then, for almost initial values Xa NxK, the following is satisfied.
(a 0) Solutions of (3) are HFD if
3 (Ax( (A)L) < 3(1(W))
and LFD otherwise.
(a < 0) Let (L) and (L) be the smallest positive and biggest negative non-zero eigenvalue of L, respectively. Solutions of (3) are ----(L))FDD if
3 (Ax fa (A-(L)) < 3 (L)). Otherwise, solutions of (3) are (1 (L)))-
Proof.
The proof follows the same reasoning as the proof for the heat equation in Theorem D.5. The difference is that we now apply Lemma D.4 instead of Lemma D.3.
Therefore, the dominating frequency is either Ax(W) la (Ar(L) or d ( (Av(L) if > and A ( fa (L)) or Xı (W) fo (A+(L) if < 0.
D.2 Frequency Dominance for Numerical Approximations of the Heat Equation
For n and E R, 0, the solution of (2) at time nhh > can approximated with an explicit Euler scheme
vec(x)(n Σ
n h² (-w & veccvx(())) ("
1=
which can be further simplified via the binomial theorem as
(ec(x)(n h) = (I h (W @ La))" vec(xo).
(1I)
Hence, it holds the representation formula
vec(x)( h)
Σ Gra -1-
hA, (W) (ALL))) v. (W) & ti (L)
In this case, the dominating frequency maximizes |1 h A, f, (AxL))). When h < ||] the product h X (W)f, fa (Ax) is guaranteed be in [-1, and
1 1 h x, (W) la (M(L)) A, (W) fa (Ar()) E [0,2].
Therefore, the dominating frequency minimizes h A, (W) fa (Ar( This is the reasoning behind the next result.
Proposition D.9. Let h E IR, h 0. Consider the fractional heat equation (2) with R. Let {x(n h be the trajectory of vectors derived by approximating (2) with an explicit Euler scheme with step size h. Suppose ||W||- Then, for almost all initial values Xo
ANN (I L)
x(n
if (A1(L) < (W),
-00000
[x(11)122
0,
otherwise.
Proof. Define
(Aav >b) = arg max {|1 hλ, fW) (A(L)) (t,,.,,,,, (1...,
r.l
By the hypothesis on h, this is equivalent (Das No) arg min (), (W) f, (Ar( (1...., [I....,N]}
r.d
32

---

Therefore, (Aa, Ao) is either (A1 AN (L)) or (Ax (1), Hence,
veo(x)(n n-000 Ca.b tha (W) @ 4/b (L)
[/vec())(n))]
||ca,b va (W) @ t/b (L) 1
If the condition Ax( (W) (A)( < (W) is satisfied,
have = 1. Then by (10), the
normalized vec(x) converges to the eigenvector of & I corresponding to the largest frequency 1 - (L) = Ax L). An application of Lemma D.2 finishes the proof.
Ax (W),, (L)) < (W) not satisfied, we have N, and the other direction follows with the same argument.
Similarly to Proposition D.9 one can prove the following results for negative fractions.
Proposition D.10. Let h R. > O. Consider the fractional heat equation (2) with a Let {x(a hi) nex be the trajectory of vectors derived by approximating the solution of (2) with explicit Euler scheme with step size h. Suppose that < || The approximated solution is (1 X-L_) if
(A+ < (x(W)) fo (L)),
and (1- x1 otherwise.
Proof. The proof follows the same reasoning as the proof of Proposition D.9 by realizing that the dominating frequencies (Aa, As) are either given by (A1(W), X+ (L)) or (Xx(W). X- (L)).
D.3 Frequency Dominance for Numerical Approximations of the Schrödinger Equation
For n and h 0, the solution of (3) at time n 0 can be approximated with an explicit Euler scheme as well. Similarly to the previous section, we can write
vec(x)(n = (I+ (W & L*))" vec(xp).
and
vec(x)(n(t) Σen, cr. (1+ihd, (W) (Ar(L)))" vi (W) X 1 (L)
The dominating frequency will be discussed in the following theorem.
Proposition D.11. Let E TR, h > O. Let (x(nh)) nEN be the trajectory of vectors derived
by approximating (3) with an explicit Euler scheme with sufficiently small step size h. Sort the eigenvalues of W in ascending order with respect to their imaginary part. Then, for almost all initial values X0
AN (I-
(I-L)
L)
8 x(nh)
f 3 (MW
<
2
[x(n)]l22
0,
otherwise.
Proof. Define
(da, No) arg max { + hA, (W) fa (M(L)) r (1....,K) K) lE fI.....N]
By definition and b, for all r and it holds
1ithA fw) fo (Ao(L)) 11 (W) (W) (AL))
(12)
Hence,
vec(x)(t) t-000
Ca.biaa (W) X th (L)
[x)e(x))t)l2
|co, (W) & tb (L) 1
We continue by determining the indices a and 6. To do so, we note that (12) is equivalent to
fa (A ( (W)) - f. (As (L)) (Ac (W)
h
fo (A (L)) 12 (W) (Ax (w))
2
33

---

for all r, L, Denote by the gap
0<E min
(fo (Ar (W)) - (1) (W)))
(r.0)((a.b)
Noting that
la (Ar (L)) AA (W)12 -6 (Ao(L)) Ao (W)I) S /// | = h [W]?
la (2r (L)) P - (L)( Aa (W))* <E
gets that (12) is satisfied for h < ||| Therefore, for sufficiently small h. the domi- nating frequencies are the ones with minimal imaginary part, i.c., either fa (A1(L) 3 (Ax (W) or f, (L)) 3 (W)) If fa (A1 (L)) 3 (W)) < fo (Av(L)) 3 (Ar(W) then and the normalized vec (x) converges to the eigenvector corresponding to the smallest frequency (L). By (10), this is also the eigenvector of I @ I L corresponding the largest frequency 1 - (L) = AN L). An application of Lemma D.2 finishes the proof.
1.
Finally, we present similar result for negative powers. Proposition D.12. Let R. h > Consider the fractional Schrödinger equation (3) with a Let (x(n)) hey be the trajectory of vectors derived by approximating the solution of (3) with an explicit Euler scheme with step size h. Suppose that is sufficiently small. Sort the eigenvalues of W in ascending order with respect their imaginary part. The appraximated solution is X4 (L))-FD if
(W) fa (++(L) < Ax(WW). (L))
and (1-A-(L))-PT otherwise.
Proof. Similar to Proposition D.11, we can prove the statement by realizing that the dominating frequencies (Nas As) in (12) are either given by (4) or (L))
E Appendix for Section 5.2
We begin this section by describing the solution of general linear matrix ODEs of the form (6) in terms of the Jordan decomposition of M. This is required when M not diagonalizable. For instance, the SNA of directed graph is not in general symmetric matrix, hence, not guaranteed to be diagonalizable. We then proceed in Appendix E.1 with the proof of Theorem 5.6.
For given matrix M CNXN, the Jordan normal form is given by
M PJP-1
where P E N is invertible matrix whose columns the generalized eigenvectors of M, and J E CNN is a block-diagonal matrix with Jordan blocks along its diagonal. Denote with ..,..., Am the eigenvalues of M and with J1,.-. Jm the corresponding Jordan blocks. Let ky be the algebraic multiplicity of the eigenvalue λ and denote with (M)) (e().....)) the generalized
eigenvectors of the Jordan block J,.
We begin by giving the following well-known result, which fully characterizes the frequencies for the solution of a linear matrix ODE.
Lemma E.I. Let M PJP-1 NNN be the Jordan normal form of M. Let x : [0, T] + IR be a solution to
x'(1) Mx((), x(0) Xo-
Then, x is given by
ki
サi- (i j)l of (M),
x(t) = Σ exp (Mi)))
j
where
Xg
Pej.
and (ef: : (l....ki) LE is the standard basis satisfying Pef (M).
34

---

Proof. By (Perko, 2001, Section 1.8), the solution can be written as
kz
exp (ML) eep (J()
Pei
(Jx
Σ
1=1 i=1
where exp (Jt) diag (fexp (J: ()))))))) and
M
t
(A
exp (J+() exp (Ar(M)
21
Since exp (JI) +i=1 exp (Jr 1), focus on single Jordan block. Fix E (l....,m,, holds
P Pexp(Jit) Σ cfef
-exxp(>(MM) (ctc) + o (to) + Φt) + tot 1 ( + to + (e...)
2
= exp((((M) ( o2 (1v) + v
C + of (2. + (M) + (M +...)
t)
o8 (M).
i=1
J=1
Bringing the direct sums together, we get
サーズ
exp (M
Σ exp (Mr(M)
ji- (M).
from which the thesis follows.
In the following, we derive formula for the solution of ODEs of the form
x((t) = Mx(t)W, x(0) = Xo for diagonal matrix W CK*k and general square matrix M CNN with Jordan normal form PJP-1, By vectorizing, we obtain the equivalent linear system
(13)
vec(x))(L) = W s M ve) vec(x)(0) Then, by properties of the Kronecker product, there holds
vec(xo).
(14)
W @ M = W 8 (PJP
(I @ (W @ J) @ p-1) (I @ () JJ(
Note that (I & P) & (I & P is not the Jordan normal form of D & M. However, we can characterize the Jordan form of W M as follows.
Lemma E.2. The Jordan decompasition of W is given by W &J = PjP-1 where a block
diagonal matrix with blocks
Wi (J)
(J)
j 3
(J)
1
my (J)
and P a diagonal matrix obtained by concatenating p,,- diag (fwjatil)
ky )n=1)
35

---

Proof. As J +i==1 Js, we can focus single Jordan block. Fix / E m.]. We have
W & J, = diag
[w) }1
w J
j=1
hence, we can focus once again on a single block. Fix j E {l..... K the Jordan decomposition of WyJJ is given by P; diag
{[ n+1
Ri
and
Jnm1
1
(J)
Ji
(J)((J) 1 (J)
To verify it, compute the (n, m) element (P.J.P.")
= - (Pr) i,k
(J:) i,k (P,'))
k,m
Since P, is diagonal matrix, the only non-null entries the diagonal; therefore, and m
(Pi n,
J1 n.m
m.m
and the only non-null entries of J, are when = or + 1, hence
P
w (J)
m 12
n,n
n,
P
(P-))
wi Wj Wj 11
n,
n,++1
a+1n-1
The thesis follows from assembling the direct sums back.
Lemma E.2 leads to the following result that fully characterizes the solution of (14) in terms of the generalized eigenvectors and eigenvalues of and W.
Proposition E.3. Consider (14) with M PJP-1 and W PJP-1 where J and P are
given in Lemma E.2. The solution of (14) is
K
Rig
vec(x)
Σ Σ exp (Ar, (W)),,
(Ars (W)) ei, 8 vo2 (M),
4=1b=1
where the coefficients L, given by
kig
vec(xo) ΣΣΣ ch (I & P)Pe: @
where {et, 12 f1..... m) iE (l...., Ky) is the standard basis satisfying Pej, vt, (M).
Proof. By Lemma E.2, the eigenvalues W&M and the corresponding eigenvectors and generalized eigenvectors are
Ar (M),, e,, @ v/l (M), (Ae, (W)) -1)))++ll ei, @ dí, (M)
for 4 E (l...., (14) given by
KJ. la E (l....,m) m. and i E ( k.. Hence, by Lemma E. the solution of
m
(x)
Σ
exp (M)Ar, (W)))
4 Σ
(i - j)! (A,, (1, (W)) (% (L)),
4-1 1==1
j=l
where the coefficients 4 are given by
m
vcc(x,)
ΣΣΣ
(I P)P(( @ (c
36

---

E.1 Proof of Theorem 5.6
In the following, we reformulate and prove Theorem 5.6.
Corollary E.4. Let G be strongly connected directed graph with SNA RNNN Consider the initial value problem in (2) with diagonal channel mixing matrix W RR and a 1. Then, for almost all initial values Xg E RNNK the solution to (2) HFD if
AK (L) < (WNx
and X1 (L) is the unique eigenvalue that minimizes the real part among all eigenvalues L. the solution LFD.
L. Otherwise,
Proof. Using the notation from Proposition E.3 and its proof, we can write the solution of the vectorized form of (2) as
K vec(x)(t) = Σ
t' (N,,
Σ 1111 1-
exp (- (W)A,, (L)
(el, @ v (L)).
f=1 j)ز
As done extensively, we separate the terms corresponding to the frequency with minimal real part. This frequency dominates the exponential converges faster than polynomials for going to infinity. Consider the case Ax( W (L))((()) < A1 (W))R(AN(L))
As X (L) is unique, the product xx(W) (L)) is the unique most negative frequency. Assume without loss of generality that Ax (W) has multiplicity one. The argument does not change for higher multiplicities as the corresponding eigenvectors orthogonal since W is diagonal. Then, Ax (W)) (L) has multiplicity one, and we calculate (t) as
Ht- (i j)!
exp (-\,, (W),, (L) 8=1
12 Σ j=l
(,))) @ v
.i = =KK, exp (.)
(ex 8 (L))
(h) 1)!
(k1 1] (Ax (W) @ v/ (L))
+ K,1
j=2
-1
+ exp (-IAk
(X)] 3 vL()
K,1
(i j)!
/=1
K
kig
Σ exp (-A, (=), (W)A,, (L)t)
(w)) 8 (V)(L)
1=111=2
1=1
j=l
(xp(-tNx(W)A((()))
(ex 3 w (L))
K.1
(ki - 1)!
Σ
(Ax (Ax(W)) (ex & v((())
+ K,1
(k1 j)m
1=2
11 1
1 x k2 (Ax())) @ v(())
+
CK11
(1 - j)!
1=1
1=1
K
- (W)A exp -
c
- k1
(i (i-)) (i1, (W)) @ 0% (I)
37

---

We can then write the normalized solution as
K11
(ek @ (L))
ki カー
Σ
(Ax(W))=> @ eL)
C11
(ki 1)!
(ki - j)!
==2
(i إ}ز (AK (W)( @ el
K
m
fi-j-ki
c1 (2 - j)! (M1, (W)) @ ((,(1)
12=2
K11
c Σ (ox & (L))
(ki j)ز (k (W)) @ v (L)
(k1 1)!
カー
إزز (An, (W)) & (L)
CK11
j=l
K
exp
(-((A, (L) -x (L)))
1=11==2
kig
ドーjート: j)! (Ar, (w)) @ v/ (L))
All summands, except the first, converge to for t going to infinity. Hence,
vt) ve)(x)(1)ll
t1+00
K.1
h
(ex
@ ((((L)
chi (ki - (ex 2 )v(())
(ki - 1)!
We apply Lemma D.2 to finish the proof for the HFD case. Note that $ (L) is an eigenvector corresponding X1 The LFD case equivalent. By Perron-Frobenius for irreducible non- negative matrices, there is no other eigenvalue with the same real part as 1 1-Ax(1 = A1 - L). O
Remark E.S. If the hypotheses are met, the convergence result also holds for L With the same reasoning, can prove that the normalized solution converges to the eigenvector corresponding to the eigenvalue of L with minimal real part. It suffices to consider the eigenvalues and generalized eigenvectors of L However, we do know the relationship between the singular values of L where we defined the fractional Laplacian, and the eigenvalues of L. Hence, it is much more challenging draw conclusions on the Dirichlet energy.
E.2 Explicit Euler
In this subsection, we show that the convergence properties of the Dirichlet energy from Theorem 5.6 are also satisfied when (2) is approximated via explicit Euler scheme.
As noted in (11), the vectorized solution (2) can be written as vec(x)(n (I h (W & L))" vec(xo),
when a = 1. We thus aim to analyze the Jordan decomposition of L" for L E C" and n E Let L= PJP-1 where the Jordan form, and P is invertible matrix of generalized eigenvectors. Consider Jordan block J, associated with the eigenvalue (MMM). For positive integer n, the n-th power of the Jordan block can be computed as:
38

---

(L)
(L)
- Ar +++1
1
(L)
2) Xr( 12++
J (L)
1
AL)
We compute the n-th power of L as L (PJP-1y) PJ and expand X
m
h
X
Pc
where iof: :E (......r... (1,.... m}} } is the standard basis and Pej (L) are the general-
ized eigenvectors of L. is easy to see that
lu
kz
L = PJ
Pej
PJ
cej
As J +i=1 JP, we can focus single Jordan block. Fix / E {I,...,,, m. and compute
PJ
Σ c/e
P (f((M)) +P ( +
(0) AM -ot + ( (LL 170 + A((M)
+
We summarize findings in the following lemma.
ki Lemma E.6. For any PJP-1 NXN and X of (L), we have
m min (k
L
Σ
n
We proceed with the main result of this subsection.
Proposition F.7. Let G be strongly connected directed graph with SNA E RNNN Consider the initial value problem in (2) with diagonal channel mixing matrix W E RKKK and a = 1. Approximate the solution (2) with an explicit Euler scheme with sufficiently small step size Then, for almost all initial values E N the following holds. If (L) is unique and
AK (L) < A1 (W)RAN (L),
(15)
the approximated solution is HFD. Otherwise, the solution is LFD.
Proof. As noted in (11), the vectorized solution to (2) with a can be written as (n (I- h (W & L))" vec(xo).
Consider the Jordan decomposition of L = PJP 1 and the Jordan decomposition of W 8 J PJP-1, where J and are specified in Lemma E.2. Then,
(nc(x)(nh) = (I + Aw @ (PJP (PJP-1))* vec(xo)
(I & P)(I hw @ J)" & P) vvec(xo)
39

---

(I& P) hPJP (1 @ P) vec(xo)
(I @ P 1 P) vec(xo)
= (I & P)P(I- - -)P(I-AJ))(( (J)" P)P) vec(xo).
Now, decompose X, into the basis of generalized eigenvectors, i.c.
K
m
voc(xo) = Σ
d( ()P)(, 8 o (I)).
1=
Then, by Lemma E.6, we have
K m
min k1y.1-1)
vec(x)(n
1=1 i i-j) (1- hA,, (W)A,,
(A,, (w)) bh (W) @ v%((L).
Now, consider the maximal frequency, i.e.,
L1, arg max (11 hA, (W)),, (L))).
41
Then, the solution vec(x) h) can be written as
min Kig 1)
i {1m)) (1 hir,, (W),
Σ Σ Σ
of G2 U1 (W) 8 v%,((L)
(1 - hdz, (Z) (L))
(16)
min ( k (-1)
h1, (L)
ΣΣ
Σ i-j {{-;)
n
(1- (W)A,, (L)) L th, (W) @ vVL(L).
With a similar argument as in the proof of Theorem 5.6, we can then see that
(ec(x)(nh)) (x)(x))n1)]
c, H2, (W) @ vL, (L) lict,,,, vr, (W) @ v1,, (L)
t1000
where wh, (L) is the eigenvector corresponding to AL. (L). Note that for almost all x, we have L,,. * O. Then w1, (L) is also an eigenvector of I L corresponding to the eigenvalue 1 Ar, (L). By Lemma D.2, we have that the approximated solution is (1 λL₃ -FD)
We finish the proof by showing that L2 if (15) satisfied, and L N otherwise. First, note that either Xx (W)RA, or A (W) (L) are the most negative real parts among all (A)AAR (ee (I....,, (I.....),, K). ry( Assume first that Ax (W))RAA has the most negative real
part, i.e., (15) holds. Then, define
max AK (W) Ar (W)))RAA
and assume h < [||] Now is easy to see that
2Λ W RA (L) hAx ()) () 22 - (L))
which is equivalent to [K, 1) = (L1, Lz). Hence, the dynamics are (1- A1 (L))--D.. As (1 -1- X (L)) highest frequency I - L, we get HFD dynamics. Similary, we can show that if X ( is the most negative frequency. we get LFD dynamics. Note that for the HFD argument, we must assume that A₁ (L) the unique eigenvalue with the smallest real part. For the L.FD argument, already given that Ay (L) has multiplicity one by Perron-Frobenius Theorem.
40

---

E.3 GCN oversmooths
Proposition E.8. Let be a strongly connected and aperiodic directed graph with SNA L R A GCN with the update rule
Xt+1 Lx, W;
where RNNK are the input node features, always oversmooths.
Proof. The proof follows similarly to the proof Proposition E.7. The difference is that instead of (16), we can write the node features after layers
K
min Aig (-1)
Σ Σ i- {--,) (A,, (W))r, (L)) (L)) or,, (W) & v/. (L).
vec(x,)
4=11=11
Now note that by Perron-Firoenius the eigenvalue Ax (L) with the largest absolute value is real and has multiplicity one. Then, maxh, (W) (L) is attained at either A1 (W)))y (L) or show that the corresponding GCN is Av (L)-FD. Now (L = A1 L), and A1 - L)-FD corresponds to LFD, hence the GCN oversmooths.
A (W)x(() (L). Equivalently to the proof of Proposition E.7,
F Appendix for the Cycle Graph Example
Consider the cycle graph with nodes numbered from 0 to - 1. Since each node has degree the SNA L A/2 is circulant matrix produced by the vector (e) eN-1) 1) /2. Denote w = exp (2mi (2r//N) NN), the eigenvectors can be computed as
1 (1,w,,e w0 ((-)))))
jj
asociated to the eigenvalue X; = cos(2Tr/N) First, we can note that X = AN-j for all j E (1...., N/2); therefore, the multiplicity of each eigenvalue is 2 except Åo and, if N is even, AN/2- Since the original matrix is symmetric, there exists a basis of real eigenvectors. A simple calculation
LAV, iLSv, = Lv, X3v; AyRV, idy
shows that RV, and Sv defined as
1
2TT n2 N
N-1 SV, VN
2
N-1
Rv,
cos
sin
VN
n=0
n=0
are two eigenvectors of the same eigenvalue Xy. To show that they are linearly independent, we compute under which conditions
0= Rv, ; + DSV,
We note that the previous condition implies that for all & N/,
(2Tjjn n
2mj n
0 coS
bsin
N
N
22
b
Va² a2 + 62 sin
arctan
N
Suppose a, b b 0, then must be
2 N
b arctan
-
+
Kn, t: E E
à
which is equivalent to
b a N, & E
arctan
2ÿ =
41

---

The left-hand side is always integer, while the right-hand hand side is an integer if and only if = 0. This reduces the conditions
2
a Cos
0
|a] sin
N
which is true if and only if 0. Consider now an even number of nodes N; the eigenspace of
AN/2 = - is
Vx/2 = VN ((-1)")>=
hence, the maximal eigenvector of I guarantees homophily 0. Consider now number of nodes N divisible by 4; the eigenspace of AN/4 = O has basis
1 VN
N-1 n=
RV
cos
SVN/4
sin
2
n=
VN
2
Their sum is then equivalent to
m
N-1
RV + SVN/4
CO8
+ sin
VN
2
2
n=
V2
m
N-1
sin
+
VN
2
n=0
2
N-1
sin
( (2n
1)
n=0
VN
hence, the mid cigenvector of guarantces homophily 1/2. A visual explanation shown Figure
42


Reviews:
[{"review_id": "1C1VK3MJyd", "reviewer": "Reviewer_ZW5p", "summary": "The authors propose a neural ODE (FLODE) that uses a fractional graph Laplacian as an alternative to GNNs, which famously suffers from oversmoothing after a few layers. The heat equation $x'(t) = -L^{\\alpha}x(t)W$ is shown to possess nice qualities wrt Dirichlet energy such that it does not always end up with a low frequency solution (i.e. oversmoothing as $t \\to \\infty$). Several experiments are done with synthetic, real, directed, and undirected graphs, and the proposed method are shown to perform comparably with existing method.", "questions": "1. If Film is a heterophilic graph, I expected the learned exponent to be closer to 0, which is the case with Squirrel and Chameleon. Do the authors have a hypothesis about why this might be the case?\n\n2. Similarly, do the authors have a hypothesis about why FLODE was not in the top three performing models for Film, Pubmed and Cora?\n\n3. Have the authors actually inspected the eigenvalues and singular values of their models and compared them to the theoretical results?\n\n4. I recommend using boldface and italics in addition to color in Table 1. In Figure 3, the line for FLODE should be thicker or have a different line style.", "limitations": "Yes.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Even though I am not an expert in neural-ODE based GNNs, I found the paper clearly written and easy to follow. The strength of this paper is in its theoretical analysis of the spectral values of their solutions. Based on a cursory reading of Di Giovanni et al (2022), and given that the learned $\\alpha$s are often not 1 (Appendix Table 4), this does not seem to be a trivial extension of existing work.", "weaknesses": "Please see questions below.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for the time taken to carefully assess our work and for the valuable feedback. We address each point individually. “W/Q” numbers the weakness or question, followed by our response.\n\n---\n> **Q1**: “If Film is a heterophilic graph, I expected the learned exponent to be closer to 0, which is the case with Squirrel and Chameleon. Do the authors have a hypothesis about why this might be the case?”\n\nFilm is a notably challenging dataset, where GNN-based models often struggle to significantly outperform models that do not incorporate the graph structure, such as MLPs. We tried to initialize the exponent $\\alpha$ to smaller values, however, the learning algorithm did not adapt $\\alpha$, and also the performance did not improve. There could be several reasons for this.\n1. Film is by far the sparsest graph among the datasets we used, which could potentially impede the propagation of information.\n2. The features in this dataset may be less expressive (with 932 dimensions compared to 2325 and 2089 for Chameleon and Squirrel, respectively).\n3. It is a class-imbalanced graph, which makes predictions more challenging.\n4. The Film dataset represents a subgraph of a larger real-world Wikipedia graph, which may result in a loss of critical information.\n\nNevertheless, these are mere hypotheses and the definitive cause remains unclear.\n\n---\n> **Q2**: “Similarly, do the authors have a hypothesis about why FLODE was not in the top three performing models for Film, Pubmed and Cora?”\n\nPlease note that FLODE is among the best models on the directed version of the Film dataset, see Tab. 2. For Pubmed and Cora, we suspect that our model's performance may have been affected by our limited computational resources and a consequently restricted hyperparameter search. Importantly, our model essentially reduces to GRAFF when we set $\\alpha=1$. Therefore, given the same hyperparameters, we should always be able to match or outperform GRAFF. We are optimistic that with more extensive tuning, our model would improve its standings on these datasets.\n\nFinally, we want to mention that Pubmed and Cora are highly homophilic datasets, where standard GNNs already perform very well as the 1-hop neighborhood is already very expressive. Hence, we do not expect improvements by using FLODE as long-range dependencies are probably not that important.\n\n---\n> **Q3**: “Have the authors actually inspected the eigenvalues and singular values of their models and compared them to the theoretical results?”\n\nWe have examined the eigenvalues of our models and compared them to our theoretical predictions, see Point 2 in our general response. These detailed results are included in the supplementary pdf, see Tab. 1 and Fig 1. The observed parameters align perfectly with our theoretical results, indicating that our method can adaptively learn to exhibit HFD behavior, particularly in heterophilic settings. This supports our claim that our proposed model dynamically adjusts its behavior according to the data.\n\n---\n> **Q4**: “I recommend using boldface and italics in addition to color in Tab. 1. In Fig. 3, the line for FLODE should be thicker or have a different line style.”\n\nWe appreciate your suggestions. We will use underline, boldface, and italics in addition to colors in Tab. 1 to ensure its accessibility. Also, in Fig. 3, we will make the line for FLODE thicker for better visibility.\n\n---\nWe value your thorough review and positive feedback on our work, especially given the depth of our theoretical analysis, which you highlighted. In light of our clarifications and the strengths you've identified, we believe that our work represents a significant advancement in addressing the oversmoothing issue in GNNs, backed by a solid theoretical foundation. We hope that our explanations and additional ablations further underscore the connection between our strong experimental results and our theoretical framework, convincing you even more."}, {"review_id": "TEVBSC2esu", "reviewer": "Reviewer_dNLj", "summary": "The authors consider the problem of classification on attributed graphs with a geometrical approach. They introduce a fractional graph Laplacian for undirected and directed graphs. They generalize the notion of Dirichlet energy to directed graphs. They study the ODEs based on this Laplacian and prove that their solutions converge to high-, low- or middle-frequency patterns, for undirected graphs and for directed graphs in a particular case. Finally they implement these ODEs as GNNs and assess their performances.", "questions": "On formatting: in the current state of the manuscript many intra-links are broken. Some figures of interest for understanding are not in the main text (I did not look at them). Maybe for table 1 the authors could use font formatting (bold, italic, ...) instead of colors.\n\nL. 245 « In accordance with the results in Section 5, we select W as a diagonal matrix. » could the authors develop?\n\nPart 6, for the directed SBM, what are the features?\n\nOut of curiosity: what happens if W is untied, i.e. taken time-dependent or different at each layer? Does one have the learned alpha small (or negative) for heterophilic datasets?", "limitations": "The authors address a main limitation of their model, its computational complexity. They may conduct an experiment showing that truncating the SVD is correct.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "I do not know well this topic and can hardly evaluate how significant are the improvements this article brings.", "weaknesses": "To my understanding a main point of this article is fractional Laplacian for directed graphs, but theorem 5.5 is restricted to the SNL.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for the time taken to carefully assess our work and for the valuable feedback. We address each point individually. “W/Q” numbers the weakness or question, followed by our response.\n\n---\n**W1**: We agree that Thm. 5.5 seemed to be restricted, and we have since made improvements to address this concern. We have now generalized Thm. 5.5 to encompass not only $\\alpha=1$, but any $\\alpha$ as long as the underlying graph has a normal SNA. This was feasible without altering our proofs since any normal SNA is unitary diagonalizable. See also Point 3 in our general response.\nMoreover, we would like to clarify some points:\n\n1. Our results not only hold for directed graphs but also for undirected. Moreover, the improvement in performance can also be seen in undirected graphs, see ablation studies or Point 2 in general response.\n\n2. The ability to capture long-range dependencies via the FGL isn't tied solely to Thm. 5.5. The FGL does also add “virtual” edges between long-distant nodes, as described in Section 4. Thus, our approach exhibits its effectiveness on tasks with (directed) heterophilic graphs. We further substantiate this claim with an additional ablation study reported in the supplementary pdf.\n\n---\n**Q1**: We apologize for the inconvenience caused by the broken intra-links. They currently refer to the appendix, and this issue will be resolved in the camera-ready version, as the main manuscript and appendix can be combined into one document.\n\nRegarding the Figs. you've mentioned, we would greatly appreciate if you could specify which ones you are referring to, so we could potentially move them to the main text, considering the availability of an extra page for the camera-ready version. Based on the suggestions of the other Reviewers, we move our ablation studies to the main paper, however, we are open for discussion if you believe other Figs. are even more important for understanding our points.  Lastly, we thank you for your suggestion about the color usage in Tab. 1. We will underline the first place, bold for the second, and italic for the third, to ensure better readability.\n\n---\n**Q2**: In Thm. 5.3 and Thm. 5.5, our theoretical findings are valid for symmetric and diagonal channel mixing matrices respectively. To align with these theoretical insights and ensure their applicability to our model, we selected $\\mathbf{W}$ to be a diagonal matrix in all our experiments. It's worth noting that all diagonal matrices inherently possess symmetric properties, which makes them suitable for satisfying the conditions of both Thms. We will clarify this in our revised paper.\n\n---\n**Q3**: For the directed SBM, we followed the original setup used in the MagNet paper, where the features are one-dimensional and are sampled from the standard Gaussian distribution.\n\nMaking $\\mathbf{W}$ time-dependent or distinct at each layer is certainly intriguing. We conducted a set of experiments with this modification to our model and found that it yielded comparable performance.  Moreover, the exponent does display similar tendencies. Specifically, for heterophilic graphs (Chameleon and Squirrel), the learned exponent is small. Meanwhile, for homophilic graphs, the learned exponent remains close to $1$. Notably, Film presents an anomaly — there's a drop in accuracy, and its learned exponent doesn't diminish, despite it being a heterophilic graph.\n\n\n| | | Film | Squirrel | Chameleon |Citeseer | Pubmed | Cora |\n|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|\n|Undirected| Accuracy |34.72±1.21|63.49±1.66|73.42±2.0|78.19±2.55|88.96±0.35|86.66±1.09|\n| | $\\alpha$ |0.92±0.01|0.30±0.04|0.38±0.06|0.97±0.007|0.98±0.05|0.95±0.02|\n|Directed|Accuracy|35.27±1.33|73.29±1.62|78.79±1.86|-|-|-|\n| | $\\alpha$ |0.94±0.004|0.25±0.08|0.25±0.08|-|-|-|\n\n---\nHowever, our theoretical analysis regarding OS would not directly extend to this adjusted model. Analyzing solutions for time-dependent ODEs is more challenging. Since one of our goals in this work was to maintain a tight relationship between theory and practice, we chose to share $\\mathbf{W}$ among all layers. Nevertheless, we are happy to include our experimental results using an untied $\\mathbf{W}$ in the appendix and identify the investigation of such time-dependent neural ODEs in both theory and practice as a promising avenue for future work. Finally, to facilitate further research in this direction, we will update our code to include an ``--no_sharing`` option. This will allow other researchers to easily explore these questions using our open-source code.\n\n---\n**Limitations**: By employing a truncated SVD and retaining $k$ singular values, we achieve the optimal $k$-rank approximation of the fractional Laplacian, as supported by the Eckart–Young–Mirsky theorem. As shown in Appendix A, Figure 8, even when using a truncated SVD that captures only 24% of the singular values, we still attain the same SOTA performance on Chameleon as with the full SVD. In the revised version of our paper, we will clarify this point and provide a reference for the corresponding theoretical result.\n\n---\nWe appreciate your thorough feedback and have incorporated all suggested changes. Given our revisions and the excellent rating you provided for the presentation, we kindly ask you to reconsider the scores for soundness, contribution, and overall score. We emphasize our novel approach in introducing a FGL, generalizing Dirichlet energy to directed graphs, ensuring that our models can mitigate OS and that our strong experimental results are fully backed up by the theoretical ones. We hope that this convinces you of our strengths."}, {"review_id": "Gzxy9iDFsP", "reviewer": "Reviewer_R3EG", "summary": "The paper extends the notion of Dirichlet energy to define oversmoothing in directed graphs. It introduces a fractional Laplacian operator that is used to create graph neural ODE architectures. The resulting designs are non-local and can alleviate oversmoothing. The paper provides empirical results (accuracy) on nine real-world benchmarks for node classification.\n", "questions": "**Questions and comments**\n\n1. Why is it called Fractional graph laplacian if it is based on the decomposition of the normalized adjacency (SNA)?\n2. Do the authors have any intuition to explain why the proposed method works so well on the Squirrel dataset?\n3. Some parts could be improved for precision and clarity, e.g.,\n    - In the paragraph 'Homophily and Heterophily', the paper uses $i \\in V$ to denote the distribution used in the expectation. What does it mean (e.g., uniform dist. over the nodes)?\n    - In the introduction, the paper says: \"a GNN that minimizes the Dirichlet energy is expected to perform well on homophilic graphs...\". However, there is no dependency on the GNN --- x is defined as initial node features.\n5. Running the ablation study on Appendix F3 to more datasets would be interesting and helpful to validate some intuitions behind the proposal.\n6. There is a vast literature on methods for tackling oversmoothing in GNNs (DropEdge, etc.), enabling such approaches for heterophilic datasets. Shouldn't these approaches also be considered as baselines for comparison?\n7. Have you tried employing an NN instead of using just a channel mixing matrix?", "limitations": "The authors discuss the limitations of the proposal, mainly mentioning the computation cost associated with SVD.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "* Overall, the paper is well-written and easy to follow. Also, the proposed framework (FLODE) is simple;\n* The proposed method works relatively well for both directed and undirected graphs (flexibility).\n", "weaknesses": "- Novelty: I rank the novelty as moderate. First, directed symmetrically normalized laplacian has been previously used. Moreover, it seems that the advantages come mainly from the choice of diffusion operator. However, it is clear that the choice of diffusion operator impacts oversmoothing, and we can design, for instance, high-pass (or non-local) graph filters to alleviate oversmoothing. Last, there is no novelty in the evaluation setup.\n- Theory: Some theoretical results focus on the properties of the fractional Laplacian, which appears only complementary to the paper. As a result, the paper does not discuss empirical results in light of the theoretical ones.\n- Experiments: The empirical results are not strong. Except for Squirrel, the performance gains (over the second-best baseline) do not seem to be statistically significant. Also, the proposed model is the best-performing model only in 4 out of 9 datasets. In addition, the paper does not consider large-scale datasets from OGB, for instance.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for the time taken to carefully assess our work and for the valuable feedback. We address each point individually. “W/Q” numbers the weakness or question, followed by our response.\n\n---\n**W1**: Please note that our work is the first in-depth theoretical analysis of the directed SNA, tying it to Dirichlet energy. We introduce a new def. for the fractional graph Laplacian (FGL) specific to directed graphs and present novel theoretical results.\n\nRegarding your comment on the choice of diffusion operator: we agree that the operator selection can have an impact on OS. However, it's not a foregone conclusion that just any diffusion operator can alleviate this issue. Through our work, we've demonstrated, both theoretically and empirically, that the FGL *combined* with the ODE framework is capable of mitigating OS. However, this does not necessarily hold for the FGL without the ODE framework. To further clarify this point, we have included new additional ablation studies (see Point 1 in our general response, Figure 1 and Table 1 in the supplementary pdf), complementing the one we conducted on the Chameleon dataset (refer to Appendix F.3).\n\n---\n**W2**: Our theoretical insights are not mere supplements but rather fundamental components of our work. They underline how the FGL, used in conjunction with graph ODEs, is adept at capturing long-range dependencies in both directed and undirected graphs. This conclusion stems from Thms. 4.1 and 5.3-5.5; the theoretical results from Section 3 are fundamental for the proofs of the above-mentioned thms.\n\nTo further emphasize this point, we have made efforts to closely associate our empirical results with our theoretical ones. See Points 1 and 2 in our general response.\n\n---\n**W3**: We value your feedback but have a different perspective on our empirical results. Our method excels on 4 out of 9 datasets and performs well especially on heterophilic graphs, ranking best on 3/6 and top three on 5/6. Additionally, our performance on the directed datasets (Fig. 3) is competitive with the SOTA.\n\nWhat distinguishes our approach is its simplicity. Unlike many models, ours doesn't rely on layer normalization, positional encodings, jumping knowledge, rewiring, dropout, non-linearity in the GNN, or extensive hyperparameter searches. Our ablation study underscores how our model's performance stems from our simple design choices.\n\nOn large-scale datasets, while our method holds potential for scalability through truncated SVD, our hardware constrains us from effectively handling datasets with 100K+ nodes. Nevertheless, our current results and theory pave the way for future advancements, especially around scaling.\n\n---\n**Q1**:  The term \"Fractional Graph Laplacian\" refers to our approach to applying the fractional power of a Laplacian. The term \"Laplacian\" is often used in graph theory to denote a matrix representation of a graph. Our approach of taking fractional powers in the singular value domain generalizes to any graph laplacian, not only the SNA or the SNL, hence the name. We hope this clarifies our choice, and we remain open to alternative suggestions!\n\n---\n**Q2**:  The Squirrel dataset is heterophilic, where nodes of different classes connect. Standard GNNs often underperform on such datasets due to OS. Our method mitigates this, explaining our enhanced performance over baselines on Squirrel. Yet, OS isn't the sole challenge in graph learning; \"oversquashing\" is another, preventing long-range node information transfer. [Song et al., 2023] suggest Squirrel experiences oversquashing. Our model's success might also stem from the virtual edges introduced by the FGL, possibly addressing oversquashing. However, this remains speculative without formal proof.\n\nAdditionally, we present ablation studies on “Squirrel”, see Point 1 in the general response.\n\n[Song et al., 2023] Song, Yunchong, et al. \"Ordered GNN: Ordering Message Passing to Deal with Heterophily and Over-smoothing.\" (2023).\n\n---\n**Q3**:  Thank you for your feedback. To address your concerns:\na) Yes, we mean a uniform distribution over nodes. We've now opted for a simple sum over the nodes to avoid confusion.\n\nb) We clarify the GNN's dependency by stating: \"A GNN whose output node features minimize the Dirichlet energy is expected to perform well on homophilic graphs.\"\n\nIf there are other areas needing clarity, please let us know. Your assistance is vital in ensuring that our work is both comprehensive and accessible to all readers.\n\n---\n**Q4**:  We appreciate your suggestion. In response, we've included a more comprehensive ablation study for the Squirrel dataset, as detailed in the supplementary pdf (see Tab. 1). This extended analysis should provide additional validation for our method.\n\nAdditionally, we added ablation studies that show that our method can scale to large depths and validate our theoretical results (see Points 1-2 in the general response and Fig. 1 in the supplementary pdf).\n\n---\n**Q5**:  See Point 3 in our general response.\n\n---\n**Q6**:  We have explored other learnable parameters beyond a diagonal channel mixing matrix, including full matrices and time-dependent matrices. These did not have the same theoretical results or empirical performances as our method, which is why we discarded them. However, we ran the experiments, see also our response to Reviewer dNLj.\n\nWhile it would be possible to apply a neural network, this approach can not be analyzed in our current framework. We believe this could be a direction for future work and we'd be glad to mention this possibility in our conclusion.\n\n---\nWe have addressed each of your comments and feel our paper stands on solid theoretical and empirical grounds. We would appreciate if you could specify the parts you consider “unsound”, given your score of 2, so we can further refine our work. If our responses have solved your concerns and you believe our paper has improved, we kindly ask you to improve the score.\n"}, {"review_id": "GDep8OBH19", "reviewer": "Reviewer_PWaE", "summary": "The authors propose two novel Fractional Graph Laplacian (FGL)-based neural Ordinary Differential Equations (ODEs): the fractional heat equation and the fractional Schrödinger equation. These solutions provide enhanced flexibility in the convergence of the Dirichlet energy and make the exponent of the fractional graph Laplacian a learnable parameter. This allows the network to adaptively decide the optimal exponent for a specific task and graph. The fractional graph Laplacian operator generalizes the Laplacian operator. The experimental results highlight the improvements of using fractional graph Laplacians, but the benefit is limited.", "questions": "1. In Table 1, could you provide the learned values for $\\lambda_{K}(W)$ and $\\lambda_{1}(W)$ for each dataset? Moreover, could you offer a detailed analysis of how these results relate to Theorem 5.3? Also, all the learned $\\alpha$ values in Table 4 are greater than 0, so according to part(i) of Theorem 5.3, the FGL-ODE is still either L-FD or H-FD. Could you clarify under which circumstances the proposed $\\lambda$-FD could be applied?\n2. The paper doesn't clearly illustrate how the fractional Laplacian operator mitigates the oversmoothing issue, both from theoretical and experimental perspectives. More in-depth explanations and experiments are needed.\n3. Is the FGL-based Neural ODE applicable to large-scale graph datasets, such as those in references [1] and [2]?\n4. The results for the Cora dataset presented in Table 1 are inferior to MLP. Could you also display the results for the homophilic dataset using the data split method from reference [5]?\n5. Table 1 seems to miss important baselines for heterophilic graphs, such as ACM-GNN[3] and ACMP[4]. It would be beneficial to include these.\n6. The paper lacks ablation studies. Could you provide the node classification results for $\\alpha=1$ in equations (4) and (5) of FLODE? As the main contribution appears to be the FGL, are there any ablation studies that use layer-wise FGL-based GNNs instead of the FGL Neural ODE?    \n7. Other questions please refer to weakness part.   \n\n  \n\n\n\n\n\n[1].Hu W, Fey M, Zitnik M, et al. Open graph benchmark: Datasets for machine learning on graphs[J]. Advances in neural information processing systems, 2020, 33: 22118-22133.  \n[2].Lim D, Hohne F, Li X, et al. Large scale learning on non-homophilous graphs: New benchmarks and strong simple methods[J]. Advances in Neural Information Processing Systems, 2021, 34: 20887-20902.  \n[3].Luan S, Hua C, Lu Q, et al. Revisiting heterophily for graph neural networks[J]. arXiv preprint arXiv:2210.07606, 2022.  \n[4].Wang Y, Yi K, Liu X, et al. ACMP: Allen-Cahn Message Passing for Graph Neural Networks with Particle Phase Transition[J]. arXiv preprint arXiv:2206.05437, 2022.  \n[5]. Shchur O, Mumme M, Bojchevski A, et al. Pitfalls of graph neural network evaluation[J]. arXiv preprint arXiv:1811.05868, 2018.", "limitations": "Yes,  the authors addressed the limitations.", "rating": 3, "confidence": 5, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The proposed FGL Neural ODE can be $\\lambda-FD$ which extends the Neural ODE-based GNNs that are limited to being either LFD or HFD.\n2. This paper generalizes the concepts related oversmoothing from undirected graphs to directed graphs.\n", "weaknesses": "1. The application of the fractional graph Laplacian operator results in the loss of the graph's sparse property, transforming it into a dense even complete graph.\n2. The paper's primary contribution appears to target oversmoothing issues, as suggested by the title. However, the evidence supporting the claim that the fractional graph Laplacian can mitigate oversmoothing in GNN is insufficient.The paper posits that the theory and experimental results presented counteract the issue of oversmoothness. However, based on [1], oversmoothness can still occur even when node features exist within an invariant space that has a rank greater than one. In this context, the paper's conclusions about oversmoothness mitigation defined using Dirichlet energy appear quite constrained. How would the proposed model fare with 256 layers stacked? How does it measure against the ODE-based methods (some are also designed to address oversmoothing), such as the recent work like GRAND[2], GraphCON [3], GRAFF [4], GRAND++[5] and GREAD [6]?\n3. While the paper introduces concepts pertaining to directed graphs, the significant Theorem 5.3 only applies to undirected graphs, and Theorem 5.5 is solely restricted to the standard Laplacian where $\\alpha=1$. Although it seems that the Fractional Graph Laplacian (FGL) may bridge the gap between homophilic and heterophilic graphs, the benefits of using FGL for directed graphs aren't clearly established.  \n4. The authors introduce several definitions, notably the Dirichlet energy for a directed graph (equation 1), but fail to provide solid justification for their choice. For instance, one could question why the energy is not defined in the following way:\n$$\\sum_{i, j=1}^N a_{i, j}\\left(\\left\\|\\frac{\\mathbf{x}_i}{\\sqrt{d_i^{\\text {in }}}}-\\frac{\\mathbf{x}_j}{\\sqrt{d_j^{\\text {in }}}}\\right\\|_2^2 +\\left\\|\\frac{\\mathbf{x}_i}{\\sqrt{d_i^{\\text {out }}}}-\\frac{\\mathbf{x}_j}{\\sqrt{d_j^{\\text {out }}}}\\right\\|_2^2\\right).$$\nThe experimental results provided do not sufficiently substantiate these definitions or the theoretical assertions made in the paper.  \n5. Very poor experimental results: While the paper focuses on directed graph Laplacian, the enhancements it introduces seem to be quite narrow, especially when considering Table 2, the only table for directed graphs. The inclusion of several definitions and theories seems to add complexity to the paper without providing a substantial contribution to the GNN community. The authors also need to include the SOTA results on those datasets.  \n6. The usage of complex numbers, as in equation (5), should be justified convincingly. While this may be commonplace in pure mathematics or physics, it needs to be shown why such a model is necessary for tackling oversmoothness in GNNs. Recent developments, such as Pytorch's support for complex numbers, do not automatically warrant their use in all contexts. If equation (5) is utilized, it only increases the model's complexity and memory consumption. Therefore, a robust justification is needed.    \n\n[1] Oono, Kenta, and Taiji Suzuki. \"Graph neural networks exponentially lose expressive power for node classification.\" arXiv preprint arXiv:1905.10947 (2019).\n\n \n\n[2] Chamberlain, Ben, et al. \"Grand: Graph neural diffusion.\" International Conference on Machine Learning. PMLR, 2021.\n\n \n\n[3] Rusch, T. Konstantin, et al. \"Graph-coupled oscillator networks.\" International Conference on Machine Learning. PMLR, 2022.\n\n \n\n[4] Di Giovanni, F., Rowbottom, J., Chamberlain, B. P., Markovich, T., and Bronstein, M. M. (2022). Graph Neural Networks as Gradient Flows: Understanding Graph Convolutions via Energy. arXiv: 2206.10991 [cs, stat].\n\n \n\n[5] Thorpe, Matthew, et al. \"GRAND++: Graph neural diffusion with a source term.\" International Conference on Learning Representations. 2021.\n\n \n\n[6] Choi, Jeongwhan, et al. \"GREAD: Graph Neural Reaction-Diffusion Networks.\" (2023).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are very grateful to the reviewer for the time taken to carefully assess our work and for the valuable feedback. We address each point individually. “W/Q” numbers the weakness or question, followed by our response.\n\n---\n\n**W1**: Using FGLs does not lead to a loss of sparsity across all graphs: the graph's density incrementally increases as $\\alpha$ goes to 0, see Appendix A, Fig. 6. Moreover, we've demonstrated in Thm. 4.1 that the FGL still respects the topology of the underlying graph. We also argue that adding some “virtual” edges can be advantageous for heterophilic graphs as nodes far apart in the input graph can directly interact.\n\n---\n\n**W2**: The notion of oversmoothing (OS) in [1] and the Dirichlet energy are intimately connected, see [Cai et al., 2020]. When the invariant subspace is the kernel of the SNL this corresponds to OS (wrt the decay of the Dirichlet energy). While the invariant subspace is always the kernel in standard GNNs, we show that this is not true for our method! This mitigates *provably* OS. To more directly address your query, Fig. 1 in the supplementary pdf illustrates that the performance of our model does not deteriorate with 256 layers. Please also note that Dirichlet energy is a standard notion to measure OS, see, e.g., [3,4] that you mention.\n\nRegarding the comparisons, [2,4] are already included in Tab. 1, and we add [3,5,6] to the already discussed models designed to mitigate OS, see also Point 4 in the general response.\n\n---\n\n**W3**: We extend Thm. 5.5, see Point 3 in our general response.\nThe efficacy of the FGL for directed graphs is substantiated in two key ways. 1) FGLs exhibit the ability to directly capture long-range dependencies by adding “virtual” edges, see also our answer to **W1**, 2) we present ablation studies on the effect of the FGL, see Point 1 in the general response. All ablation studies show that the combination of the FGL and the ODE framework improves the performance on directed graphs.\n\nWe hope that this is enough evidence for the benefits of using FGL in directed graphs.\n\n---\n\n**W4**: The justification for our def. of Dirichlet energy is actually elaborated on lines 103 ff., immediately after its introduction. Additionally, our defs. of the Dirichlet energy (Def. 3.1), SNA, and SNL (Def. 3.2) and their relationship to the Dirichlet energy (Prop 3.6) in directed graphs recovers the standard defs. in the undirected case, see, e.g., [Cai et al., 2020].\n\n---\n\n**W5**: We would like to respectfully clarify a few points. In addition to Tab. 2, Fig. 3 also presents results for 11 directed graphs.\nOur model demonstrates superior performance on 4 out of 9 real-world datasets, and is on pair with the SOTA on the 11 directed graphs in Fig. 3 (see Tab. 5 in Appendix F.2). Moreover, our approach does not rely heavily on hyperparam. tuning and avoids unnecessary complexity in the GNN layers (jumping knowledge, non-linearities, dropout, or positional encodings), maintaining the comprehensibility and relevance of our theoretical results. This is in contrast to the models we compare with, e.g., GREAD, GraphCON, GRAND.\n\nWe believe our contribution to the GNN community lies in our novel approach to understanding directed graphs and fractional Laplacians, an area that has been less explored in the past. Our work can provide a foundation for further exploration and advancements in this area.\n\nRegarding SOTA results, we welcome suggestions on methods to include in our revision.\n\n---\n\n**W6**: Our justification is straightforward: complex numbers are natural when dealing with directed graphs, as eigenvalues and eigenvectors can be complex. Also one the ODEs we consider, the Schrödinger equation, is inherently complex. Moreover, the complexity of the GNN layers scales linearly to the hidden dimension, i.e., complex numbers are not problematic at all. Thus, the use of complex numbers in our work is rooted in necessity and applicability, rather than the availability of software support. This was already used by related works on directed graphs, e.g., MagNet [Zhang et al., 2021].\n\n---\n\n**Q1**: Please see Point 2 in the general response for the first two questions. Regarding $\\lambda$-FD: as noted in Fig. 2, graphs with medium homophily could benefit from this. One could for example initialze the exponent with a negative value when the homophily level is known in prior. Our theory offers the backbone for future work to analyze these use-cases.\n\n---\n\n**Q2**: Please see our response to **W2** and Point 2 in the general response. Also, note that we actually *prove* that our approach mitigates OS, see Section 5.\n\n---\n\n**Q3**: Our approach can be scaled using methods such as randomized or truncated SVD. We noted this in Fig. 8  in Appendix A, where 20% of the singular values are already sufficient for SOTA performance.\nHowever,  given our present resource constraints we are currently unable to run experiments on such large-scale datasets. We acknowledge this as a significant future direction for our work, though it is beyond the immediate scope of this study.\n\n---\n\n**Q4**: Cora is relatively simple as both, the nodes’ features and their 1-hop neighbourhood, are expressive. However, as you rightly suggested, employing fewer training nodes (as in [5]) does indeed showcase the utility of our method, even without any hyperparameter tuning:\n\n| model | Cora | Citeseer | Pubmed |\n|:-:|:-:|:-:|:-:|\n| MLP |  58.2±2.1  |  59.1±2.3  |  70.0±2.1  |\n| Ours | 77.8±1.1 | 69.1±1.6 | 74.2±2.2 |\n\nWe happily incorporate these findings into the paper's appendix.\n\n---\n\n**Q5**: Please see Point 4 in the general response.\n\n---\n\n**Q6**: Please see Point 1 in the general response.\n\n---\n\nWe've thoroughly addressed your concerns. Considering the positive feedback from other reviewers and our comprehensive revisions, we respectfully request a reevaluation of your score. Specifically, if any sections still appear unsound or lacking in presentation, please guide us to them.\n"}]

Author rebuttal:
We thank the reviewers for their thorough and insightful remarks. We fully implemented all remarks in the revised version of the paper, which we believe improved the paper significantly.

**1. More extensive ablation studies**

To augment our ablation studies on Chameleon and Citeseer (see Appendix F.3), we perform several new experiments, the results of which can be found in the new supplementary pdf. First, we carried out a comprehensive ablation study on “Squirrel”, see Tab. 2. Here, we systematically remove various parts of our model (the learnable exponent, the ODE, both, and the complex weight matrix). The results show a significant drop in performance when any of these components are removed, reinforcing the fact that our model's convincing performance is indeed a result of the combination of the fractional Laplacian and the neural ODE framework.

Upon request of the Reviewers, we conduct another ablation study to investigate the role of depth on Chameleon, Citeseer, Cora, and Squirrel datasets. The results are depicted in Fig. 1 of the supplementary pdf. The results demonstrate that the neural ODE framework enables GNNs to scale to large depths (256 layers). Finally, we see that the fractional Laplacian improves over the standard Laplacian in the heterophilic graphs which is supported by our claims in the main paper. We highlight that using only the fractional Laplacian without the neural ODE framework oftentimes outperforms the standard Laplacian with the neural ODE framework, see ablations in Tab. 2 in the supplementary pdf and Tab. 7 in Appendix F.3. This indicates the importance of the long-range connections built by the fractional Laplacian.

Please also be aware that Cora, Citeseer, and Pubmed are graphs with high homophily levels. While we do not anticipate improved performance on these graphs (as capturing long-range dependencies is less important on the graphs), we also do not expect a performance drop. For instance, we could simply set $\alpha=1$, which should ensure that we achieve at least the same performance as GRAFF, provided the hyperparameters are equal.

**2. Highlighting the intimate connection between our theoretical and experimental findings**

We further demonstrate the close alignment of our theoretical and experimental results in Fig. 1 of the supplementary pdf, which enables us to precisely anticipate when the models will exhibit HFD or LFD behaviors. In this context, we calculated parameters (according to Thm. 5.3) and illustrated at each depth the expected and observed behaviors. For Squirrel and Chameleon, which are heterophilic graphs, we observe that both their theoretical and empirical behaviors are HFD. Additionally, the learned exponent is small. In contrast, for Cora and Citeseer, we see the opposite.

We have also incorporated the calculated parameters (according to Thm. 5.3) for Tab. 1 in the main (see Tab. 1 in the supplementary pdf). Here, we employ the best hyperparams to solve both fractional heat and Schrödinger graph ODEs., further substantiating the intimate link between our theoretical advancements and practical applications. This tight connection provides strong validation for our theoretical framework, indicating its utility in effectively guiding and predicting the behavior of GNNs in practice.

Please also note that we could enforce LFD/HFD behavior by using diagonal channel mixing matrices with only positive/negative entries. However, this doesn’t lead to the same empirical results which may be due to the limited expressivity of the resulting network.

**3. Generalization of Thm. 5.3**

We appreciate the insights provided by two of the reviewers and generalize Thm. 5.3 to any directed graph with normal SNA. The updated version of the Thm. will be included in the revised paper, and the existing proof remains valid as any normal matrix is unitarily diagonalizable, allowing the same proof technique used in Thm. 5.3. We regard further generalizations to arbitrary directed graphs as an intriguing open question for future exploration.

We thank the reviewers for inspiring this straightforward but crucial extension of Thm. 5.3.

**4. Inclusion of baseline models for experiments**

In response to the insightful suggestions from several reviewers, we have expanded our set of compared models (see Tab. below) to include recent neural ODE-based GNNs and other methods designed to mitigate oversmoothing (OS). These additions encompass GREAD, GraphCon, ACMP, and GCN/GAT with DropEdge. Please note that for models not evaluated using the same (standard) splits, we have run their code on these standard splits, a detail we will highlight with a footnote in our revised paper.

The results will be added in our paper:

| Model | Film | Squirrel | Chameleon | Citeseer | Pubmed | Cora |
| - | - | - | - | - | - | - |
| GREAD [Choi et al., 2023] | **37.90±1.17** | 59.22±1.44 | 71.38±1.3 | 77.60±1.81 | **90.23±0.55** | **88.57±0.66** |
| GraphCon [Rusch et al., 2022] | 35.58±1.24 | 35.51±1.40 | 49.63±1.89 | 76.36±2.67 | 88.01±0.47 | 87.22±1.48 |
| ACMP [Wang et al., 2023] | 34.93±1.26 | 40.05±1.53 | 57.59±2.09 | 76.71±1.77 | 87.79±0.47 | 87.71±0.95 |
| GCN+DropEdge [Rong et al., 2020] | 29.93±0.80 | 41.30±1.77 | 59.06±2.04 | 76.57±2.68 | 86.97±0.42 | 83.54±1.06 |
| GAT+DropEdge [Rong et al., 2020] | 28.95 ± 0.76 | 41.27±1.76 | 58.95±2.13 | 76.13±2.20 | 86.91±0.45 | 83.54±1.06 |
| Ours (This work) | 37.16±1.42 | **64.23±1.84** | **73.60±1.55** | **78.07±1.62** | 89.02±0.38 | 86.44±1.17 |

We also move the paragraph detailing our baseline models from Appendix F.1 (lines 704 ff.) to the main text, improving clarity on our selection of robust baseline models for heterophilic graphs. Given the reviewers' emphasis on baseline models addressing oversmoothing — a concern already addressed in our initial submission — this shift seemed important.


---
We address the remaining concerns as part of the response to each Reviewer.


Meta-review: Most reviewers recommended the paper to be accepted.  Reviewer PWaE engaged in a detailed and lengthy discussion with the reviewers and maintained firmly that the paper should be rejected. 
I found the discussion between the authors and the reviewers interesting, and the authors' rebuttal sufficient for the acceptance of the paper. Nevertheless, I strongly recommend that the authors add the additional experiments on larger datasets as reviewer PWaE suggested.
Decision: Accept (poster)