Graph Neural Networks as Gradient Flows: understanding graph convolutions via energy
Authors: Francesco Di Giovanni
James Rowbottom
Benjamin Paul Chamberlain
Thomas Markovich
Michael M. Bronstein
Venue: Submitted to ICLR 2023 (2023)
Keywords: Graph Neural Networks
Gradient flows
Energy functionals
Spectral theory

Abstract: Gradient flows are differential equations that minimize an energy functional and constitute the main descriptors of physical systems. We apply this formalism to Graph Neural Networks (GNNs) to develop new frameworks for learning on graphs as well as provide a better theoretical understanding of existing ones. We derive GNNs as a gradient flow equation of a parametric energy that provides a physics-inspired interpretation of GNNs as learning particle dynamics in the feature space. In particular, we show that in graph convolutional models (GCN), the positive/negative eigenvalues of the channel mixing matrix correspond to attractive/repulsive forces between adjacent features. We rigorously prove how the channel-mixing can learn to steer the dynamics towards low or high frequencies, which allows to deal with heterophilic graphs. We show that the same class of energies is decreasing along a larger family of GNNs; albeit not gradient flows, they retain their inductive bias. We experimentally evaluate an instance of the gradient flow framework that is principled, more efficient than GCN, and achieves competitive performance on graph datasets of varying homophily often outperforming recent baselines specifically designed to target heterophily.

Under review as a conference paper at ICLR 2023
GRAPH NEURAL NETWORKS AS GRADIENT FLOWS: UNDERSTANDING GRAPH CONVOLUTIONS VIA ENERGY
Anonymous authors Paper under double-blind review
ABSTRACT
Gradient flows are differential equations that minimize an energy functional and constitute the main descriptors of physical systems. We apply this formalism to Graph Neural Networks (GNNs) to develop new frameworks for learning on graphs as well as provide better theoretical understanding of existing ones. We derive GNNs as a gradient flow equation of parametric energy that provides physics-inspired interpretation of GNNs learning particle dynamics in the feature space. In particular, we show that in graph convolutional models (GCN), the positive/negative eigenvalues of the channel mixing matrix correspond tractive/repulsive forces between adjacent features. We rigorously prove how the channel-mixing can learn to steer the dynamics towards low high frequencies, which allows to deal with heterophilic graphs. We show that the same class of energies is decreasing along a larger family of GNNs; albeit not gradient flows, they retain their inductive bias. We experimentally evaluate instance of the gradient flow framework that is principled, more efficient than GCN, and achieves com- petitive performance graph datasets of varying homophily often outperforming recent baselines specifically designed to target heterophily.
INTRODUCTION
Graph neural networks (GNNs) Sperduti 1993
Goller & Kuchler 1996 Gori et al. 2005 Scarselli
al. 2008 Bruna
2014 Defferrard al. 2016 Kipf & Welling 201 7 Battaglia lia al. 2016
Gilmer et al. 2017 7) have become the standard ML tool for dealing with different types of relations and interactions. Limitations of GNNs that have recently attracted attention in the literature are over-smotthing (node features becoming increasingly similar with the depth of the model, see Nt & Maehara (2019): [Oono Suzuki (2020): Cai & Wang (2020): [Zhou et al. (202T), over-squashing (the difficulty of message passing propagate information the graph, see Alon & Yahav (202 I); Topping et | (2022) ), and poor performance on heterophilic data (i.c. where adjacent nodes tend to
have different labeles, see Pei al. 2020); Zhu et al. (2020); Bo et al. (021): Yan et al. (2021). )
General motivations and contributions. In the spirit of neural ODEs (Haber & Ruthotto 2018 Chen al.] 2018), regard (residual) GNNs as discrete dynamical systems. A fundamental idea in physics is that particles evolve by minimizing an energy: one can then study the dynamics through the functional expression of the energy. The class of differential equations that minimize an energy are called gradient flows and their extension and analysis in the context of GNNs represent the main focus of this work. We study two ways of understanding the dynamics induced by GNNs: starting from the energy functional or from the evolution equations.
time
From energy to evolution equations: a new conceptual approach to GNNs. We propose a general framework where one parame- terises energy functional and then takes the GNN equations follow the direction of steepest descent such energy. We introduce class of energy functionals that extend those adopted for label
Figure 1: Gradient flow dynam- ics: attractive and repulsive forces lead to a process able to separate heterophilic labels.

---

Under review as a conference paper at ICLR 2023
propagation (Zhou & Schölkopf 2005) and whose gradient flow equations consist of generalized graph convolutions (GCN -ype architectures, (Kipf & Welling 2017 with symmetric weights. We provide a physical interpretation for GNNs as multi-particle dynamics: this new framework sheds light on the role of the 'channel-mixing' matrix used in graph convolutional models as an edge-wise potential inducing attraction (repulsion) via its positive (negative) eigenvalues. We conduct theoretical analysis of the dynamics including explicit expansions of the GNN learned features, showing that differently from other continuous models, the gradient flow can learn to magnify either the low or high frequencies. This also establishes new links to techniques like residual connections and negative edge weights that have been previously used in heterophilic settings. We experimentally evaluate our framework using gradient flow equations yielding principled variant of GCNs that is also more efficient due to weight symmetry and sharing across layers. Our experiments demonstrate competitive performance on homophilic and heterophilic graphs of varying size.
From evolution equations to energy: understanding graph convolutions via multi-particle dynamics. Recent works of Cai & Wang (2020); Bodnar al. (2022) studied the behaviour of the Dirichlet energy in graph convolutional models in order to determine if (over)smoothing of the features is occurring. The key idea is that the monotonicity of an energy functional along a system of equations conveys significant information about the dynamics, both in terms of its dominating effects and limit points. However, these results are restricted to the classical Dirichlet energy and assume (non-residual) graph convolutions activated by the ReLU nonlinearity. We extend this approach by proving that a much more general multi-particon convolutions with symmetric weights and with respect functions, Our result sheds light onto the dynamics of non-linear graph convolutions showing that the 'channel-mixe mixxng" matrix used in GCN-type models can be interpreted as potential in feature space that promotes alignment (repulsion) of adjacent node features depending on its spectrum.
energy is in fact decreasing along residual graph-
more general class of nonlinear activation
Outline. In Section we review non- con-paraettric instances of gradient flows on graphs: the heat equation and label propagation. In Section B we extend this approach to the parametric case by introducing class of energies that generalize the one used for label propagation and whose associated gradient flows are continuous graph convolutions. We provide physical interpretation for the energy showing that it can induce attraction and repulsion along edges. In Section 4 we discretize the gradient flow into GNN update equations and derive explicit expansions of the learned node representations highlighting how the spectrum of the channel-mixing W controls whether the dynamics is dominated by the low high frequencies of the graph Laplacian. To knowledge, ours is the first analysis that studies the interplay of the spectral properties of the graph Laplacian and the channel mixing matrix. In Section 5 we extend the theory by showing that the same multi-particle energy introduced in Section still decreases along more general graph convolutions with symmetric weights, meaning that the physics interpretation is preserved. In Section 6 we evaluate the framework for node classification on a broad range of datasets.
Related work. Our analysis is related to studying GNNs as filters (Defferrard et al. 2016 Hammond
ct al. 2019 Balcilar et al.] 2020 He et al. 2021) and adopts techniques similar to Oono & Suzuki
(2020): Cai & Wang (2220) Gradient flows were adapted from geometry (Eells & Sampson 1964),
to image processing (Kimmel et al. 1997). label propagation (Zhou & Schölkopf 2005) and recently
in ML (Sander et al. 2022) for the analysis of Transformers (Vaswani et al. 2017 7. Our work follows
the spirit of GNNs as continuous dynamical systems (Xhonneux al. 2020 Zang & Wang 2020
Chamberlain et al. 202 Eliasof et al. 2021 Chamberlain et al. 202 Bodnar et al. 2022
Rusch
al. 2022).
Notations. Let G EV be an undirected graph with n nodes. Its adjacency matrix A is defined as d.j = (f) E E and zero otherwise. We let D = diag(d (d)) be the degree matrix and define the normalized adjacency A == D ADD/AAD-1 11/2, We denote by F E Rnxd the matrix of d-dimensional node features, by f. E Rd its i-th row (transposed), by f" IR" its r-th column, and by vec(F) E Bdd the vectorization of F obtained by stacking its columns. Given a symmetric matrix B, we let XB AB
denote its most positive and negative eigenvalues, respectively, and PB be its spectral radius. f(t) denotes the temporal derivative, & is the Kronecker product and 'a.e.' means almost every w.r.t Lebesgue measure. Proofs and additional results appear in the Appendix.
2

---

Under review as a conference paper at ICLR 2023
2 GRADIENT FLOWS ON GRAPHS: THE NON-PARAMETRIC CASE
In this Section we review important concepts on graphs and two examples of non-parametric gradient flows that partly motivate our approach to the GNN framework.
What is gradient flow? Consider an N --dimensional dynamical system governed by the evolution
equation F(t) = (DE(F(t) that evolves some initial state F(0) for time t 2 0. In deep learning,
the discretisation of such differential equations using the Euler method allows draw a parallel
between iterations of numerical solver and the layers of neural network (Haber & Ruthotto 2018)
Chen al., 2018). We say that the evolution equation is a gradient flow if there exists an energy functional : TR + R such that ODE(F()) = (E(F(()). Since (F) = (F
the energy & is decreasing along the solution F(t) of such equations. The existence of & and the knowledge of its functional expression allow for a better understanding of the dynamical system.
A prototypical gradient flow: heat equation. Let F Rnxd be the matrix representation of vector features assigned to each node in G. Its graph gradient is defined edge-wise as (VF); :=
f/Vdd - f t/Vdd We can then set the Laplacian as A := div v (the divergence div is the adjoint of V), represented by A = I- À - 0. We refer to the eigenvalues of A as frequencies: the lowest frequency is always while the highest frequency Pa S 2 (Chung & Graham 1997). The heat equation on each channel is the system f"
Af" (t), for 1 S S d. This an example of gradient flow: if we stack the columns of F into vec(F) E Rnd, we can rewrite the heat equation as
rec(F((t) Ve (vec(((()()))
where Dir IRnd nd
+ R is the (graph) Dirichlet energy defined by Zhou & Schölkopf 2005
8 Dir (F)
1
Σ ((VF)bl (IJ)EE
trace( F (F) (Fec(F), (Ia & A) Veec
4
Dir measures the smoothness of the signal since it accounts for the variations of F (i.e., its gradient) on the edges. The evolution of F()) by the heat equation I1 decreases the Dirichlet energy E (F(t)) in the limit Dir (F + o)) = 0, attained by the projection of the initial state F(0) onto ker A)
more general gradient flow: label propagation. Assume we have a graph G, featuressures Fo and labels {y:} on Virsin C V, and that we want predict the labels on Veest V. Zhou & Schölkopi (2005) proposed label propagation (LP) where they first extend the input labels (0) outside the training set as y, (0) for each i \ Vtrain and then solve the equation:
(1) ( - ( (O)).
This is another example of gradient flow; in fact. Zhou & Schölkopi (2005) originally introduced the following energy and then derived the aforementioned update formula l in order to minimize it:
(()) = VELP (Y))), ELP (Y) = Dir (Y) + 4 Y(O)
(3)
The prediction is then attained by the signal that minimizes both Dir which enforces smoothness and the fitting term arising from the available labels (a form of soft boundary conditions).
Motivations and goals. In graph ML problems, we often also have node features that can be leveraged for the label prediction. Our goal is to extend the gradient flow formalism from the porammraettric case (heat equation and label propagation) to deep learning setting, where we (i) parameterise an energy functional and let the GNN equations be the associated gradient flow, and (ii) investigate when existing GNNs admit an energy that is decreasing along their evolution equations.
GRADIENT FLOWS ON GRAPHS: THE PARAMETRIC CASE
We can think of (residual) graph neural network a parametric evolution equation F(t) GNNo F()) discretized using the Fuler method with fixed time step 0 < < < 1:
F(L + r) = F()) + GGNNe F()).
(4)
3

---

Under review as a conference paper at ICLR 2023
Each iteration corresponds to GNN layer, which in general can have different set of parameters (tt). We choose GNN, to be the gradient flow of some parametric class of energies Eg R - R generalizing Dir resulting in feature evolution by F(t) = - VE (t()) starting from input features F(O), with {0} learned via backpropagation on the task Ioss function. This approach extends the LP technique to a framework where the parameters we learn can be interpreted as "finding the right notion of smoothness' for our task. In fact, minimizing Et. as in Equation (3) works only the labels are smooth an assumption known as homophily. We investigate how learning more general energy yields gradient flow GNNs that can also perform well on heterophilic data.
3.1 ENERGIES GIVING RISE TO GRAPI-CONVOLUTIONNL MODELS
Similarly to the LP approach in Equation (3 our first step consists in choosing parametric class of energy functionals {Ee) giving rise to the GNN equations via gradient flow. GNNs of the convolutional flavor (Bronstein et al. 202T) evolve the features via (4) using some parametric rule GNN,(G, Fo) typically consisting of two operations: applying a shared linear transformation to the features ('channel mixing') and propagating them along the edges ('diffusion') "). Accordingly, introduce the class of (generalized) graph convolutions:
(5) where the learnable parameters (0()) are the d x d weight matrices I2 Wt, and W, acting on each node feature vector independently and performing channel mixing: the normalized adjacency A performs the diffusion of features from adjacent nodes. The setting T = 1, no residual connection, and Ω2 W; 0 corresponds to GCN (Kipf& & Welling 2017)). The case of S2 + 0 results in an anisotropic instance of GraphS AGE Hamilton et al. (2017), while by choosing 22 = 0 and W and
7) = F(t) + (-F(t) + ( F(O)W.)
-FO)).
W, as convex combinations with the identity we recover GCNII (Chen et al. 2020
We consider a class of energies {Ee} consisting of quadratic terms
Eg (F)
2 Σ (f, St.))
A (fs, Wi,) + (F
(F. F(O)),
(6)
2
Eoot parameterised by d x d weight matrices Ω W. We motivate our choice by first recovering the
ghair
non-parametric cases of Section Ω = w = I, and sp" -00=00 0, then Eg = gDir as per Equation (2)
choosing 8d as Le-penalty gives &g ELP per Equation (3). We can also recover manifold
harmonic energies applied to graphs (see Appendix B.5 More importantly, if (F, (())
Σ, (f,, (O)),, for W E Rdddd we can rewrite
(F) (vec(F), J & La
W A)vec(F) (W I,)vec (F(O))
(7)
and then derive its gradient flow as:
Ω 5 S
W + WT 2
Ftt) = VPE (F(()) = F(()
+ AF(()
F(O)W.
(8)
2
Since S2 W appear in Equation (8) in a symmetrized way, without loss of generality we can assume Ω and W to be symmetric d x d channel mixing matrices. Therefore, Equation ( simplifies
Ftt) = F()) + AF - FOO)W.
Thus, quadratic energy as Equation (35) leads to continuous linear graph convolutions with symmetric weights shared over time. Equivalently, for generalized graph convolutions to fit the gradient flow formalism, the channel-mixing matrices must be symmetric. Importantly, while reducing the number of parameters and offering gradient flow interpretation of the GNN, this symmetric constraint does not diminish its expressive power (Hu et al. 2019). Next, we show that E has simple interpretation in terms of pairwise forces among adjacent features.
3.2 ATTRACTION AND REPULSSON: PHYSICS- FRAMEWORK
Why gradient flows? A multi-particle point of view. Consider the node features as particles in d with energy Eg. The first term gext is independent of the pairwise interactions and hence represents an 'external' energy in the feature space. The second term Epair instead accounts for pairwise interactions along edges via the symmetric matrix W and hence represents 'internal'
4

---

Under review as a conference paper at ICLR 2023
energy. We set the source term 40 to zero and write W = o e_e_, _, by decomposing it into components with positive and negative eigenvalues. We can then rewrite &g in Equation (6) as
1 Σ (C. (f. ()C) + EII Eo (F) = 2
( (10)
attraction which we have derived in Appendix To understand the dynamics induced by the minimization of Eo by the gradient flow (36), recall that the edge gradient (VF);; measures the difference between features f; and f,. We note that: (i) if S commutes with W, then the projections of (VF),; onto ker () remain invariant and are preserved along the gradient flow; (ii) The channel-mixing W encodes attractive edge-wise interactions via its positive-definite component ++ since the gradient terms ||0 (VF)) decrease along the solution of Equation (36), hence resulting in W
dampening
repulsion
smoothing effect where adjacent node features f; and f, are 'aligned'; (ii) The channel-mixing encodes repulsive edge-wise interactions via its negative-definite component O_ since the gradient terms ||@ (VF);s increase along the solution of Equation (36). hence resulting in a sharpening effect which could be desirable on heterophilic graphs where we need to disentangle adjacent node representations. Next, we formalize the smoothing vs sharpening effects by introducing a new quantity to monitor along GNN to assess whether the latter is magnifying the low or high frequencies.
Low vs high frequency enhancement. Attractive forces minimize the edge gradients and are associated with smoothing effects which magnify low frequencies, while repulsive forces increase the edge gradients and hence afford sharpening action enhancing the high frequencies. Since we are interested in finding which frequency is dominating the dynamics, we monitor the Dirichlet energy along the normalized solution: E (() This is the Rayleigh quotient of La @ A and so it satisfies 0 I. EDr (F)/1 S pa/2 (see Appendix A.2). If the normalized Dirichlet energy is approaching its minimum, then the lowest frequency component is dominating, whereas if the normalized Dirichlet energy is converging to its maximum, then the dynamics is dominated by the highest frequencies. This allows to introduce the following
Definition 3.1. F(t) GNN ece (G, F(t)) initialized F(0) is Low/High-Frequencyy Dommnnnt + P2/2) for t ++ co.
(L/HFD) if Dir (F(t)/||F((II
0 (respectively, g (F())
we provide justifications and explicit examples. If a graph is homophilic, we expect a smoothing or LFD dynamics enhancing the low-frequency components to be successful for node classification (Wu et al., 2019: Klicpera et al.] (019). In the opposite case of heterophily, the high-frequency components might contain more relevant information for separating classes (Bo et al., 2021) - the classical example being the eigenvector of A with largest frequency e separating bipartite graph. Accordingly. an ideal framework for learning on graphs must at least accommodate both of these opposite scenarios by being able to induce either an LFD or HFD dynamics.
In Appendix B.
We can now investigate the gradient flow equations of the energy in Equation (LO).
Theorem 3.2 (Informal). The continuous gradient flow in Equation (36) can learn to be either LFD (mostly edge-wise altractive) or HFD (mostly edge-wise repulsive) depending on the spectrum of W.
(A precise result, along with convergence rates, is stated Theorem B.3 in the Appendix). Informally, Theorem [.2oremm.2]]soowe shows that the gradient flow in Equation (36) is expressive enough to induce repulsion along edges needed- - as expected based on the decomposition of Eg in Equation (IO). As argued above, a dynamical system that can never be HFD might instead struggle on heterophilic graphs where the feature signal needs to be sharpened rather than smoothed out. The property that the gradient flow can be HFD is con-triviia and in fact some continuous-ting GNNs - such as those introduced in Xhonneux et al. (2020); Chamberlain et al. (222 Eliasof et al. (202I) - are never HFD and as a result suffer on heterophilic datasets (as confirmed in our experiments in Table
Theorem 3.3 (Informal). Models like CGNN, GRAND and PDE- - GCN p are never HFD.
We refer to Theorem B.4 for statement including convergence rates and over-smoothing results.
Message of Section 3
We introduced an energy Eg allowing to learn attractive/repulsive forces along inducing an LFD/HFD dynamics as per Theorem 3.2 We argue that energies rather than evolution equations should be the object to parameterise for deriving more principled GNNs that are easier to interpret and analyse. This is studied next.
edges via the spectrum of the
5

---

Under review as a conference paper at ICLR 2023
4 FROM ENERGY TO EVOLUTION EQUATIONS: GNNS AS GRADIENT FLOWS
In order to connect our theory to practice, we discretize Equation 36 as in Equation replacing continuous time by fixed steps corresponding to GNN layers.
4.1 DISCRETE GRADIENT FLOWS AND SPECTRAL ANALYSIS
we use the Euler scheme with step size T to solve Equation (36. In our framework we parameterise the energy rather than the equations, which leads to symmetric thannel-mixing matrices S2, W E _dxdd The use of the explicit Euler discretization yields a residual architecture:
As in Equation
F(t+ F(t) + (((())+ A - F(O))) F(0) d'EN (Fu),
(11)
where an encoder V'EN : IR=XP + Ra processes input features Fo and the prediction ehD EDE (T))) is produced by a decoder VDDE 1n + == Here, k is the number of label classes, 7 = mt is the integration time, and m is the number of layers. We note that (i) non-linear activations can be included in VEN, dDDE making the entire model non-linear; (ii) since the framework is residual, even if the message-passing is linear, this is not equivalent to collapsing the dynamics into a single layer with diffusion matrix A m as done in Wu et al. 2019 see Equation (37 in the Appendix.
Interaction between the graph and channel-mixing spectra. We restrict our theoretical analysis to the gradient flows in Equation ( where remove dampening and source term effects (i.e., 8
W 0, which corresponds to residual GCN). Our technique consists in vectorizing the solution
F(L) + vec(F(()) and rewriting the update as ve +L)) = vec(F(L)) + (W @ A) ve)(F(())
(see Appendix A.2 for details). In particular, once we choose bases {Φ))) } and {Φf² eigenvectors for W and A respectively, we can write the solution after m layers explicitly:
} of orthonormal
n-1
voc(F(mT)) = ΣE (1+ xA)) m G.
OF
(12)
where (0) == ((O), Φ, & ΦΦ²). We see that the interaction of the spectra (AW) and {A²} the 'driving" factor for the dynamics, with positive (negative) eigenvalues of W magnifying the frequencies X < 1 (> 1 respectively). In the following we let X denote the most positive/negative eigenvalue of W with associated eigenvectors Φ W Note that Φ is the Laplacian eigenvector associated with largest frequency مم. We now consider the following:
X (ps - 1)) < | < 2(T(2 100---24 PPa))-
(13)
The first inequality means that the negative eigenvalues of W dominate the positive ones (once we factor in the graph spectrum contribution), while the second a constraint on the step-size since i is too large, then we no longer approximate the gradient flow in Equation (36).
Theorem 4.1. Let m be the number of layers. Consider F(t + 7) = F()) + T AF()), with
symmetric W. If Equation (3) holds, then there exists 6 < s.t. for all i E V we have:
((()) (1+ (pp 1)) (c-m-1(0) @ - (i) Φw )0
(14)
Conversely, if XW (pa 1)) (1 then
(((r)) (1 + > m (c+.(()) Vd. Φh +8())
(15)
We report the explicit value of 8 in Equation (38) in Appendix C.1 We now comment on the consequences of Theorem 4.1 Equation (4) implies that if the negative eigenvalues of W are sufficiently larger than the positive ones (in absolute value, as per Equation (1h), then repulsive forces and hence high frequencies dominate. Indeed for é E we have f (mT) Φ { &w at the fastest scale, up to lower order terms in the number of layers. Thus feature representation (((T becomes dominated by a multiple of &w E Rdd only depending on the value taken by the Laplacian eigenvector Φ at node i. On the other hand, if Equation (55) holds, then at the largest scale we have f, (mer) ~ Vdd Φw E Rdd, meaning that the node representation becomes dominated by a multiple of &f only depending on the degree of i- which recovers the over-smoothing phenomenon (Nr & Maehara 2019 Oono Suzuki 2020).
(1)((()))) increase the depth, any
OOr arguments extend trivially to the degenerate case.

---

Under review as a conference paper at ICLR 2023
Corollary 4.2. If Equation
14 holds, then the system is HFD for a.e. FFFT) + Fxc s.I.
F(O) and
sc PaPx for each r. Conversely, if Equation 15) holds, then the system is LFD for a.e. (() and F(mT)/ [ Foc s.t. A = 0 for each
Remark. In general neither the highest nor the lowest frequency Laplacian eigenvectors constitute ideal classifiers and in fact we always have finite depth so that F F(m7) also depends on the lower- order terms of the asymptotic expansion in Theorem 4.1 Whether the dynamics is LFD or HFD will affect if the lower or higher frequencies have a larger contribution to the prediction; indeed, we can compute the 'impact' of each graph frequency explicitly thanks to Equation (12).
4.2 CONNECTIONS TO EXISTING RESULTS
Residual connection. The following result shows that the residual connection is crucial:
Theorem 4.3. If not bipartite, and we remove the residual connection, i.e. F(t+r) F. with W symmetric, then the dynamics is LFD for a.e. F()) independent of the spectrum of W.
Differently from previous over-smoothing results of Oono & Suzuki (2020); Cai & Wang 2020),
here we have no constraints on the spectral radius of W coming from the graph topology. In other words, the residual connection fully enables the channel-mixing to steer the evolution towards low or high frequencies depending on the task. If we drop the residual connection, W is less powerful; this is also confirmed by our ablation studies (see Figure B in the Appendix).
W A W ( be the eigendecomposition of W yielding the Fourier coefficients ((t ) = F(I)&W W We rewrite the discretized gradient flow F( T) (tt) 7 AF W in the Fourier domain of W as Z(t r) = Z(t) + 7 (t) W and note that along the eigenvectors of W. if AW < O then the dynamics is equivalent to flipping the sign of the edges. This shows that negative edge weight mechanisms proposed Li al. (2020): Bo et al.)
Negative eigenvalues flip the edge signs. Let W
22021); Yan et al. (2021) for heterophilic graphs can be achieved with simple GCN model where the channel-mixine matrix W has negative eigenvalues. We refer Equation (39 in the appendix for thorough discussion and derivation.
The message of Section 4 Discrete gradient flows of &g are equivalent to linear graph convolutions with symmetric weights shared across layers. This provides 'multi-particle" interpretation for graph convolutions and sheds light onto the dynamics they generate. We can derive simple expansions of the learned features and show that the interaction between the eigenvectors and spectra of W and A is what drives the dynamics and determines its dominating effects. Convolutional GNN models can deal with heterophily if the channel mixing matrix has negative eigenvalues.
FROM EVOLUTION EQUATIONS TO ENERGY: INTERPRETING GNNS VIA En
Analysing energies along GNNs is one approach to investigate their dynamics. Cai & Wang 2020
Bodnar et al. (2022) showed that Dir is decreasing (exponentially) along some classes of graph convolutions, implying mversmoothing see also Rusch et al. (2022). ). In this Section, we start from time-continuous graph convolutions as in Equation (5), with a acting elementwise:
Ftt o ( F(t) AF - F
(16)
Although this is no longer necessarily a gradient flow due to a,, we prove that if the weights are symmetric, then Eo in Equation (35) still decreases along Equation (16).
Theorem 5.1. Consider a : IR + R satisfying + zo(xx) 2 O. If solves Equation (6) with S W symmetric, then t + (F (t)) is decreasing. If we discretize the system using the Euler method step size T and C+ denotes the most positive eigenvalue of S2 & I, W & A, then
method with
(F + -)) - (F(I)) S C+ CI - F())
An important consequence of Theorem 5. is that for non-linear graph convolutions with symmetric weights, the physical interpretation is preserved since the same multi-pann particle energy Eg in Equa- tion (00) is decreasing along the solution. Note that in the discrete can be interpreted as a Lipschitz regularity result. Namely, the channel channel-mixing W still induces attraction/repulsion along edges via its positive/negative eigenvalues (more explicitly, Lemma
the energy monotonicity
7

---

Under review as a conference paper at ICLR 2023
We again emphasize that the requirement of symmetric weights is not restrictive thanks to the universal approximation results Hu et al. (219)). Theorem [5.1 differs from Cai & Wang (2020); Bodnar et al. (2022) in two ways: (i) It asserts monotomicity of energy En more general than & gr since it is parametric and in fact also able to enhance the high frequencies; and (ii) it holds for an infinite class of non-linear activations (beyond ReLU).. So far we have considered time-indee indeeeeddnnnn energies. We can generalize our discussion to energies of the form Ea(- t) whose potentials now vary in time. The equations then take the form of Equation 5 with 52, and W, symmetric.
D.I.
The message of Section 5 We show that graph convolutions with symmetric weights identify curves along which the multi particie energy Eg decreases hence acting as 'approximate" gradient flows. Despite the non-linear activation of a we can still interpret the learming dynamics of convolution on graphs as finding the 'right' edge-wise attractiie/repplliiv potentials through channel mixing.
EXPERIMENTS
In Theorem 4.1 we have shown that the subset of linear graph convolutions in Equation (5) with shared, symmetric weights is characterized by the strong inductive bias that the multi-particle energy Eg in Equation (10) being minimized along the equations. Although not gradient flow, even when we activate the equations with as in Theorem 5.1 we can preserve such inductive bias since the same energy is decreasing. This means that both frameworks can provably induce attraction or repulsion along edges thanks to the spectrum of the maxine..mixing. We validate our theoretical analysis by testing that these principled (and more efficient) classes of convolutional models along which Es decreases can compete with baselines designed to target heterophilic graphs.
The model and the parameterisation. In the following we evaluate a subclass of gradient flows in Equation [I] giving rise framework termed GRAFF (Gradient Flow Framework):
GRAFF F(t +T++T) F(t) +
( (-F(t)diag(w) W :F((0),
(17)
where w E Rdd and B E T. W is symmetric x d-matrix shared across layers and (node-wise) encoder and decoder are MLPs. We consider two possible implementations for W: (i) diagonally- dominant (see Appendix E) where we learn an off-diagonal symmetric matrix and the diagonal terms separately, and (ii) the case with W diagonal and report best numbers over these two configurations. We note that both these parameterisations allow the model to control the spectrum of W more easily which we know be essential from Theorem 4.1 we refer to the methodology description in Appendix E for further details. By Theorem 5.1 if we 'activate' Equation I17 as
GRAFFNL: : F(t + Ftt) + Ta (-F(t)diag(w) + (F(t) - (0))
(18)
with a s.t. (x)) 2 O, then Eg in Equation 35) is decreasing, so that we can think of such equations as more general "approximate" gradient flows termed GRAFFNL (where NL stands for non- non-rin
Complexity. GRAFF scales O([Vpd + m|E(d), where and d are input feature and hidden dimension respectively, with p 2 d usually, and m the number of layers. Note that GCN has complexity (m))E((p d)) and in fact our model is slightly faster than GCN mainly due to the preliminary encoding performed node-wise rather than edge-wise; main baselines on heterophilic graphs like GGCN and Sheaf learn edge-wise weights based on features which is slower (as confirmed in Figure in Appendix E). Moreover, for GRAFF the number of parameters scales as O(pd + d2) while for other baselines they scale with the number of layers least O(pd + mdf).
Synthetic experiments. To investigate our claims we first use the synthetic Cora dataset of (Zhu al.) 2020 Appendix G) where graphs are generated for target levels of homophily see Appendix EX Figure reports the test accuracy vs true label homophily. For Neg-prod we set W - W, WT so to only have non-positive eigenvalues: we see that this is better than the opposite case of prod (W- W. on low-homophily (and viceversa on high-homophilly. This confirms Theorem 4.1 where have shown that the gradient flow can be HFD that is generally desirable with fow-homophily through the negative eigenvalues of W. In practice 'non-signed" variants
gnel
nog
Teat
Torgethomapnly
Figure 2: Synthetic experiments with controlled homophily.
8

---

Under review as a conference paper at ICLR 2023
Wisconsin 0.21
Cornell
Squirrel 0.22
Chameleon 0.23
Citeseer
Pulmed
Hom lovel
0.11
0.30
0.22
0.74
0.80
0.51
#Nodes
183
25]
183
7,600
5,201
2,277
3,327
18.712
2:708
*Edges
295
466
250
26,752
198,493
31,421
4,676
44,327
5,278
*Classes
7
GGCN
84.86
4.65
86
85.08
6.63
37.54
1.56
66-17
1.58
71.14
77.14
80. 15
0.37
87.95
1.06
GPRGNN
4.36
82
80.27
34.63
31.61
1.24
46.58
77.
57.54
0.38
87.90
1.18
H2GCN
81.86
87 1.88
52. 70
28
55.70
1.00
36.48
1.86
60:11 1: 2:15
77.1
1.57
89 00-40-01
0.38
87.87
1:20
OCNI
77:57
3.83
50.30
40
77. 96
3.79
57.44 + 1.30
38.47 + 1.58
63.86
3.04
77.33
1.48
50:15
0.43
88.57 + 1.26
Gcom-GCN
00.76
2.77
64.51 3.06
001.5
+ 3.67
31.59
+ 1.15
38.15
+ 0.312
60:00
2.81
02.0022:
1.15
80.90
+ 0.47
80.35
+ 1.07
PairNorm
60:27 1: 4.34
48.43
6.14
58
3:1:15
27:40
1: 1.24
50.44
1: 2.04
62.74 1. 2.82
73.59
1.47
87 53
0.44
85.79
1: 1.01
ACII
61.18
6.01
0.00
65.45
0.50
1.04
GCN
05.14-6.18
51
00.54
5.30
27.32
1.10
53.43
64.82+222 2:244
76.50
0.50
1.27
GAT
32.16.1.6.63
49.41
1.09
61.89
5.05
27.44
0.38
10.72 1. 1.35
60:26 2.50
76.55
1.25
57.30 1.10
56.53
0.48
MLP
80.81 4.75
85.20
3:31
61.80-+6-00
96.53
+ 0.70
28.77 + 1.56
46.21
2.00
74.02
1.90
73.00
2.00
57.
0.37
CONN
71.35 4.06
74.31
7.26
06.222
7.00
35.95
+ 0.98
29.24 + 1.409
1.66 46.89
76.91
1.81
87. 20
0.40
87. 10
1.35
GRAND
75.68 1.2....225 7.25
79.
3.61
7.09 52.16
55.02
1.01
00.05
1.50
51.67
1.
76.
1.77
89.02
0.51
S7.36
0.96
Sheal (max)
15.05:t:55
1.15
1.13
GRAFF
84.05
6.30
37.11 +
58.72
0.84
71.08
1.75
77.3
50.04
0.44
1.03
GRAFFAL FaL
36.49 d: 4.84
87.20 2.52
77.30
35.96 41 0.95
59.01 4: 1.31
71.38 4: 1.47
70.81
1.12
89.81 2: 0.50
57.81 1.13
Table 1: Node-classificttion results. Top three models are coloured by First, Second, Third
like GRAFF are more flexible and outperform GCN with low homophily, confirming Theorem
where we have shown that without a residual connection convolutional models are LFD irrespectively of the spectrum of W further results in Figure ii Appendix [E]
Real world experiments. In Table I we test GRAFF and GRAFFNI NL on datasets with varying
homophily (Sen al.] 2008 RRozemerrckii et al. 2021 Pei et al. 2020) (details in Appendix [E.4). We use results provided in (Yan et al. [021] Table 1), which include GCNs models, GAT (Veličković et al. 2018), PairNorm (Zhao & Akoglu 2019) and models designed for heterophily (GGCN (Yan et al.
2021b, Geom-GCN (Pei et al. 2020). H2GCN (Zhu et al. (020) and GPRGNN (Chien et al. 2021))
For Sheaf (Bodnar et al. (022). , a recent strong baseline with heterophily, we took the best performing variant (out of six) for each dataset. We include continuous baselines CGNN (Xhonneux al. 2020)
and GRAND (Chamberlain et al. 2021a) to corroborate Theorem 3.3 Training, validation and test splits are taken from Pei et al. (2020) for all datasets for comparison. We also evaluate GRAFF on larger heterophilic datasets discussed in Lim et al. (2021) (see Appendix (E66) and we compare with
further recent baselines in Appendix E.7
Results. GRAFF and GRAFFNI. are both versions of graph convolutions with stronger 'inductive bias' given by the energy En decreasing along the solution; in fact, we can recover them from graph convolutions by simply requiring that the channel-mixing is symmetric and shared across layers. Nonetheless they achieve competitive results on all datasets often outperforming slower and more complex models. They are extremely competitive on more homophilic datasets as well, in contrast with the performance of models like Sheaf mainly designed to handle heterophily.
CONCLUSIONS
We argued that when studying and developing GNNs we should focus on energy functionals rather than the evolution equations. We introduced a new framework for GNNs where the evolution is a gradient flow of a multi-particle learnable energy. This gives rise to principled graph convolutions where the channel-mixing is symmetric matrix and induces attraction (repulsion) along edges via its positive (negative) eigenvalues. We explored the thcoretical implications by investigating the dominating terms in the learned feature expansion and corroborated that this graph convolutional framework can perform well in heterophilic settings. We proved that existing (generalized) graph convolutions maintain the dynamics induced by the same class of multi-particle energies if the channel- mixing is symmetric even when they are not strictly gradient flows due to non-linear activations; this provides deeper connection between energy functionals and GNNs and extends several recent results that have monitored the classical Dirichlet energy along GCNs to shed light on their dynamics.
L Limitations and future works. We limited our attention to a class of energy functionals whose gradient flows give rise to evolution equations of the generalized graph convolution type. In future work, we plan to study other families of energies that generalize different GNN architectures and provide new models that are more 'physics' -inspired. We will also investigate time-dependent energy functionals and how to generalize our results to this setting. To the best of our knowledge, our analysis is first step into studying the interaction of the graph and "channel-mixing spectra. In future work, we will explore other more general dynamics (i.e., that are neither LFD nor HFD).
9

---

Under review as a conference paper at ICLR 2023
REFERENCES
Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practical implications. In International Conference on Learning Representations, 2021. URL /ttpttts///peereieeee net / fform?iii 10-uun/1c-1800PhhC// 100OP
Kristen M Altenburger and Johan Ugander. Monophily in social networks introduces similarity among friends-of-friends. Nature human behaviour, 2(4):284-29) 2900, 2018.
Muhammet Balcilar, Guillaume Renton, Pierre Héroux, Benoit Gaüzère, Sébastien Adam, and Paul Honeine. Analyzing the expressive power of graph neural networks in spectral perspective. In International Conference on Learning Representations. 2020.
Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al. Interaction networks for learning about objects, relations and physics. Advances in neural information processing systems, 29,209, 2016.
Lukas Biewald. Experiment tracking with weights and biases, 2020. URL https: www. wandb com/ Software
Coftwareevaaill available from wandb.com.
Deyu Bo, Xiao Wang. Chuan Shi, and Huawei Shen. Beyond low-frequency information in graph convolutional networks. In AAAI AAAI Press, 2021.
Cristian Bodnar, Francesco Di Giovanni, Benjamin Paul Chamberlain, Pietro Liò, and Michael M Bronstein. Neural sheaf diffusion: A topological perspective on heterophily and oversmoothing in gnns. arXiv preprint wrXivv2202.04577 2022.
Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint 104.::2004413475 2021.
Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and locally connected networks on graphs. 2nd International Conference on Learning Representations, ICLR 2014, 2014.
Chen Cai and Yusu Wang. A note over-smoothing for graph neural networks. arXiv preprint arXiv:2006. 13318, 2020.
Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. Grand: Graph neural diffusion. In International Conference on Machine Learning, pp. 1407-1418. PMLR, 2021a.
Benjamin Chamberlain, James Rowbottom, Davide Eynard, Francesco Di Giovanni, Xiaowen Dong, and Michael Bronstein. Beltrami flow and neural diffusion on graphs. Advances in Neural Information Processing Systems, 34, 2021b.
Ming Chen, Zhewei Wei, Zengfeng Huang, Bolin Ding, and Yaliang Li. Simple and deep graph convolutional networks. In International Conference on Machine Learning, pp. 1725-1735. PMLR, 2020.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances in neural information processing systems, 31, 2018.
Eli Chien, Jianhao Peng, Pan Li, and Olgica Milenkovic. Adaptive universal generalized pagerank graph neural network. 9th International Conference on Learning Representations, ICLR 2021, 2021.
Fan RK Chung and Fan Chung Graham. Spectral graph theory. Number 92, American Mathematical Soc., 1997.
Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. Advances in neural information processing systems, 29, 2016.
James Eells and Joseph H Sampson. Harmonic mappings of riemannian manifolds. American journaaa of mathematics, 86(1) 86(1):109-))) 1609 1964.
10

---

Under review as a conference paper at ICLR 2023
Moshe Eliasof, Eldad Haber, and Eran Treister. Pde-gen: Novel architectures for graph neural networks motivated by partial differential equations. Advances in Neural Information Processing Systems, 34, 2021.
Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR Workshop Representation Learning on Graphs and Manifolds, 2019.
Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In International Conference Machine Learning, pp- 1263-1272. PMLR. 2017.
Christoph Goller and Andreas Kuchler. Learning task-dependent distributed representations by backpropagation through structure. In Proceedings of International Conference on Neural Networks (1CNN'96), volume
pp. 347-352. IEEE, 1996.
Marco Gori, Gabriele Monfardini, and Franco Scarselli. A new model for learning in graph domains. In Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005., volume pp. 729-734. 734. IEEE, 2005.
Haber and L.. Ruthotto. Stable architectures for deep neural networks. Inverse Problems, 34, 2018.
Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. Advances in neural information processing systems, 30, 2017.
David K Hammond, Pierre Vandergheynst, and Rémi Gribonval. The spectral graph wavelet transform: Fundamental theory and fast computation. In Vertex-Frequency Analysis of Graph Signals, pp. 141-175. Springer, 2019.
Mingguo He, Zhewei Wei, Hongteng Xu, et al. Bernnet: Learning arbitrary graph spectral filters via bernstein approximation. Advances in Neural Information Processing Systems, 34, 2021.
Shell Xu Hu, Sergey Zagoruyko, and Nikos Komodakis. Exploring weight symmetry deep neural networks. Computer Vision and Image Understanding, 187:102786, 2019.
Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec. Open graph benchmark: Datasets for machine learning graphs, 2021.
Ron Kimmel, Nir Sochen, and Ravi Malladi. From high energy physics to low level vision. In International Conference Scale-Space Theories Computer Vision, pp. 236-247. Springer, 1997.
Thomas N. Kipf and Max Welling. Semi-Supervised Classification with Graph Convolutional Networks. In Proceedings of the 5th International Conference on Learning Representations, ICLR '17, 2017. URL / //pp//eeeeeeeeeeeeeeeeeeeee net/ htwparrreeieeeeeeeeoornreeeeeeee htte://ooommnndsAaaaaa
Johannes Klicpera, Stefan Weissenberger, and Stephan Günnemann. Diffusion improves graph learning. Proceedings of the 33rd International Conference on Neural Information Processing Systems, 2019.
Jure Leskovec and Andrej Krevl. Snap datasets: Stanford large network dataset collection, 2014.
Jure Leskovec, Jon Kleinberg, and Christos Faloutsos. Graphs over time: densification laws, shrinking diameters and possible explanations. In Proceedings of the eleventh ACM SIGKDD international conference on Knowledge discovery in data mining, pp. 177-187, 2005.
Yu Li, Yuan Tian, Jiawei Zhang, and Yi Chang. Learning signed network embedding via graph attention. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp. 4772-4779, 2020.
Derek Lim, Felix Hohne, Xiuyu Li, Sijia Linda Huang. Vaishnavi Gupta, Omkar Bhalerao, and Ser Nam Lim. Large scale learning on non-howe non-homophilous graphs: New benchmarks and strong simple methods. Advances in Neural Information Processing Systems, 34:20887-20902, 2021.

---

Under review as a conference paper at ICLR 2023
Vijay Lingam, Rahul Ragesh, Arun Iyer, and Sundararajan Sellamanickam. Simple truncated svd based model for node classification on heterophilic graphs. arXiv preprint TrXi:2210.11280,, 2021.
Sitao Luan, Chenqing Hua, Qincheng Lu, Jiaqi Zhu, Mingde Zhao, Shuyuan Zhang, Xiao-Wen Chang, and Doina Precup. Is heterophily a real nightmare for graph neural networks to do node classification? arXiv preprint orXiv:2109.05644 2021.
Sunil Kumar Maurya, Xin Liu, and Tsuyoshi Murata. Improving graph neural networks with simple architecture design. arXiv preprint arXiv:2 2021.
Hoang Nt and Takanori Maehara. Revisiting graph neural networks: All we have is low-pass filters. arXiv preprint wee:::9950099550 2019.
Kenta Oono and Taiji Suzuki. Graph neural networks exponentially lose expressive power for node classification. In International Conference on Learning Representations, 2020.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vito,, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang. Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. NeurlPS. 2019.
Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. Geom-gen: Geometric graph convolutional networks. In 8th International Conference on Learning Representations, ICLR 2020, 2020.
Pietro Perona and Jitendra Malik. Scale-space and edge detection using anisotropic diffusion. PAMI, 12(7))629-639) 1990.
Benedek Rozemberczki, Carl Allen, and Rik Sarkar. Multi-scale attributed node embedding. Journal of Complex Networks, 9(2)cnnab014, 2021.
T Konstantin Rusch, Benjamin P Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael M Bronstein. Graph-coupled oscillator networks. In International Conference on Machine Learning, 2022.
Michael E Sander, Pierre Ablin, Mathieu Blondel, and Gabriel Peyré. Sinkformers: Transformers with doubly stochastic attention. In International Conference on Artificial Intelligence and Statistics, pp. 3515-3530, PMLR, 2022.
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE transactions on neural networks, 20( 20(1)x61-80) 2008.
Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Ras Rad.
Collective classification network data. Al magazine, 29(3):93-93, 2008
Oleksandr Shchur, Maximilian Mumme, Aleksandar Bojchevski, and Stephan Günnemann. Pitfalls of graph neural network evaluation. In NIPS workshop. 2018.
Alessandro Sperduti. Encoding labeled graphs by labeling raam. Advances in Neural Information Processing Systems, 6, 1993.
Jake Topping. Francesco Di Giovanni, Benjamin Paul Chamberlain, Xiaowen Dong, and Michael M Bronstein. Understanding over-squashing and bottlenecks on graphs via curvature. International Conference on Learning Representations, 2022.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Eukasz Kaiser, and Illia Polosukhin. Attention all you need. Advances in neural information processing systems, 30, 2017.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In International Conference on Learning Representations, 2018. URL htll/httts://openreett http:eewwwww.eeeeeeeeeeeeeeeeeee
net.
12

---

Under review as a conference paper at ICLR 2023
Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger. Sim- plifying graph convolutional networks. In International conference on machine learning. pp. 6861-6871. PMLR, 2019.
.ml www.srpppessseeeeeee v97 wul 9e html
Louis-Pascal Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. In In-
lernational Conference on Machine Learning, pp. 10432-10441. PMLR, 2020. URL http:
//proceedings mw.. press/ /t/p//ll111////// hww11//hhoneee2.a..mm html
Yujun Yan, Milad Hashemi, Kevin Swersky, Yaoqing Yang, and Danai Koutra. Two sides of the same coin: Heterophily and oversmoothing in graph convolutional neural networks. arXiv preprint orXiv:2102.06462, 2021.
Chengxi Zang and Fei Wang. Neural dynamics on complex networks. In Proceedings of the 26th ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pp. 892-99999999 902, 2020.
Lingxiao Zhao and Leman Akoglu. Pairnorm: Tackling oversmoothing in gnns. arXiv preprint arXiv:1909.12220 2019.
Elena Zheleva and Lise Getoor. join or not to join: the illusion of privacy in social networks with mixed public and private user profiles. In Proceedings of the 18th international conference World wide web, pp. 531-540, 2009.
Dengyong Zhou and Bernhard Schölkopf. Regularization on discrete spaces. In Joint Pattern Recognition Symposium, pp. 361-368. Springer, 2005.
Kaixiong Zhou, Xiao Huang, Daochen Zha, Rui Chen, Li Soo-Hyun Choi, and Xia Hu. Dirichlet energy constrained learning for deep graph neural networks. Advances in Neural Information Processing Systems, 34:21834-21846, 2021.
Jiong Zhu, Yujun Yan, Lingxiao Zhao, Mark Heimann, Leman Akoglu, and Danai Koutra. Beyond homophily in graph neural networks: Current limitations and effective designs. Advances in Neural Information Processing Systems, 33:7793-7804, 2020.
13

---

Under review as a conference paper at ICLR 2023
OVERVIEW OF THE APPENDIX
To facilitate navigating the appendix, where we report several additional theoretical results, analysis of different cases along with further experiments and ablation studies, we provide the following detailed outline.
In Appendix A.1 we review properties of the classical Dirichlet energy on manifolds that inspired traditional PDE variational methods for image processing whose extension to graphs and GNNs more specifically partly constitutes one of the main motivations of work. We also review important elementary properties of the Kronecker product of matrices that are used throughout our proof's in Appendix |A.2 We also comment on the choice of the normalization (and symmetrization) of the graph Laplacian, briefly mentioning the impact of different choices.
In Appendix B.I] we derive the energy decomposition reported in Equation (TO). In Ap- pendix B.2 we derive additional rigorous results to justify our characterization of LFD and HFD dynamics in Definition 3.1 along with explicit examples. We also formalize more explicitly and quantitatively Theorem 3.2 in Theorem B.3 In Appendix [B.3 we report more explicit statement with convergence rates and wver-smoothing results which covers the informal version in Theorem 3.3 In Appendix [B.4 we explore the special case of Ω W which is equivalent to choosing rather than A as message-pasee arguments as to why propagating messages using A rather than the graph Laplacian is actually 'more robust'. Finally in Appendix [B.5 we formally derive an analogy between the continuous energy used for manifolds (images) and subset of the parametric energies in Equation (6)
matrix providing new
In Appendix C we prove the main results of Section namely Theorem and Theorem
Corollary
Appendix we prove Theorem 5.1 and an extra result confirming that even in the non- linear case the channel-mixing W still induces attraction and repulsion via its spectrum hence magnifying the low or high frequencies respectively.
In Appendix we report additional details on hyperparameter tuning, datasets adopted, further synthetic and ablation studies, along with extra experiments on larger heterophilic datasets in Appendix E.6
Additional notations and conventions used throughout the appendix. Any graph is taken to be connected. We order the eigenvalues of the graph Laplacian 0 A < x .... A
PA < 2 with associated orthonormal basis of eigenvectors {Ø )a-1 so that in particular AΦd O. Moreover, given symmetric matrix B, we generally denote the spectrum of B by spec(B) () and if B > 0, then gap(B) denotes the positive smallest eigenvalue of B. Finally, if we we always take the norm be the Frobenius one and tacitly assume that the dynamics is s.t. the solution is not zero.
have
write F(
A PROOFS AND ADDITIONAL DETAILS OF SECTION 2
A.1 DISCUSSION ON CONTINUOUS DIRICHLET ENERGY AND HARMONIC MAPS
Starting point: a geometric parallelism. To motivate a gradient-flow floww approach for GNNs, we start from the continuous case (see Appendix A.1 for details). Consider smooth map f IR. - (Rd h)) with h a constant metric represented by H = 0. The Dirichlet energy of f is defined by
d
( E(J, h) = gM
Σ q.r=1 1=1
|| dx
(19)
and measures the 'smoothness' of f. A natural approach to find minimizers of E called harmonic maps was introduced in Bells & Sampson (1964) and consists in studying the gradient flow of E, wherein a given map f(0) = fo is evolved according to f(t) = (f)). These type of evolution equations have historically been the core of variational and PDE-based image processing;
14

---

Under review as a conference paper at ICLR 2023
in particular, gradient flows of the Dirichlet energy were shown Kimmel al. (1997) to recover the Perona-Malin nonlinear diffusion Perona & Malik (990)..
In this subsection briefly expand the formulation of continuous Dirichlet energy in Section to provide more context. Consider a smooth map f (M. - h), where N is usually larger manifold embed M into, and g, h are Riemannian metrics domain and codomain respectively. The Dirichlet energy of f is defined by
E(Sg, g. h)
with ldflg the norm of the Jacobian of f measured with respect to and (M.g) is standard Fuclidean space R", N IRd and h is constant positive semi-definite matrix, then we can rewrite the Dirichlet energy in a more familiar form as
8(
trace (DS du
Σ
ΣE
, (z))d
The Dirichlet energy measures the smoothness of the map S, and indeed if h is the identity in Rd. then we recover the classical definition
E(f) = Σ (I
Gradient flow of Dirichlet energy. Minimizers of & referred to harmonic maps important objects in geometry: to mention a few, geodesics, minimal isometric immersions and maps S M Fd solving 4,f = 0 are all instances of harmonic maps. To identify such critical points, one computes the first variation of the energy E along an arbitrary direction 0z which can be written as
(0.f) = =- (T.(f),
for some tensor field with explicit form
(Tyae (f)) == Agn fa + hy 1 0.f Of
for 1 S a S dim(N), with {yn) local coordinates on N and r Christoffel symbols. It follows that harmonic maps identified by the condition (f)) 0. In Eells & Sampson (1964), the pivotal idea of harmonic map flow which has shaped much of modern research in geometric analysis was
introduced for the first time: in order identify minimizers of an input map fo is evolved along the direction of (minus) the gradient of the energy & leading to the dynamics
a, = (f)).
(20)
As a special case, when the target space the classical Euclidean space one recovers the heat equation induced by the input Riemannian structure. We also note that when (M, g) is a surface representing
an image and f (u1, U2) + (U1, u2 Φ(u1, u2)) with á a color map, then Equation (20) becomes
0 = div(C, V
(21)
with Cg a constant depending on the metric on M. If now let g to depend on d, one can recover the celebrated Perona-Malis flow Kimmel et al. (1997).
A.2 REVIEW OF KRONECKER PRODUCT AND PROPERTIES OF LAPLACIAN KERNEL
Kronecker product. In this subsection we summarize a few relevant notions pertaining the Kro- necker product of matrices that are going be applied throughout our spectral analysis of gradient flow equations for GNNs in both the continuous and discrete time setting.
Given matricial equation of the form
Y = AXB,
15

---

Under review as a conference paper at ICLR 2023
vectorize X and Y by stacking columns into vec(X) and vec( Y) respectively, and rewrite the previous system as
vec(YY B x A) (X)
(22)
If A and B are symmetric with spectra spec(A) A) and spec () respectively, then the spectrum of B & A is given by vectors, then A is an eigenvalue of B @ A with eigenvector xx:
Namely, Ax x and By XBy, for and y non-zero
[ 2 A) & (A 1 y xx.
(23)
One can also define the Kronecker sum of matrices A R and B E Raxd as
A + B: = A 8 L B,
(24)
with spectrum spec(A A + B) = (^^ + AB A E spec(A), A E spec
Additional details Dir and the choice of Laplacian. energy Dir is defined by
We recall that the classical graph Dirichlet
1
Dir (F)
2 Irace (FT)F),
where the (unusual) extra factor of is to avoid rescaling the gradient flow by 2 - which is the more common convention. We can use the Kronecker product to rewrite the Dirichlet energy
Dir (F)
(ec(F) (Ia @ Avvc
(25))
from which we immediately derive that Vee E and hence recover the gradient flow in Equation leading to the graph heat equation across each channel.
(F) = (I. 2 A) A)vec(F) since A symmetric
Before we further comment the characterizations of LFD and HFD dynamics, review the main choices graph Laplacian and the associated harmonic signals (i.e. how we can characterize the kemnel spaces of the given Laplacian operator). Recall that throughout the appendix always assume I A is
that the underlying graph is connected. The symmetrically normalized Laplacian symmetric, positive semi-definite with harmonic space of the form Chung & Graham (1997
ker( A := span (D+ 1,, 1, = 1)T....)))...
(26)
This confirms that if given GNN evolution F(t) (F()) with initial condition F(O)
over-smooths meaning that ( + 0 for \ - oo for each column 1 S r S d, then the only information persisting in the asymptotic regime the degree and any dependence
the input features is lost, as studied in Oono & Suzuki (2020); Cai & Wang (2020). A slightly different behaviour occurs if instead of A, we consider the unnormalized Laplacian
D A with kernel span(1,.) meaning that if Lf + 0 as t + oc for each S S d. then any node would be embedded to single point, hence making any separation task impossible. The same consequence applies to the random walk Laplacian Anw = I D- A. In particular, we note that generally a row-stochastic matrix not symmetric - it was, then this would in fact be doubly-stochastie - and the same applies to the random-walk Laplacian (a special exception is given by the class of regular graphs). In fact, in general any dynamical system governed by A Ra (or simply energy due to the lack of symmetry, as further confirmed below in Equation
A) is not the gradient flow of an
B PROOFS AND ADDITIONAL DETAILS OF SECTION
B.1 ATTRACTION VS REPULSSION; PHYSICS-INSPIRES FRAMEWORK
We first note that the system in Equation 36 can be written using the Kronecker product as
(ec(F(tt) - @ I, (el))vec(FFF) + (W & A)vec(F((t) (W & vev (Frcc(((O)) (F))..
If this is the gradient flow of F Eg (F), then we would have
EFFF 12 I, W W A,
(27)
16

---

Under review as a conference paper at ICLR 2023
which must be symmetric due to the Hessian of a function being symmetric. The latter means
(S2 (2) In (WT @
which is satisfied if and only if both Q2 and W symmetric. This shows that Equation (3) is the gradient flow of Eo if and only if Qt and W are symmetric.
We now rely the spectral decomposition of W to rewrite & explicitly in terms of attractive and
repulsive interactions. If we have spectral decomposition W = &w W (BW)² we can separate the positive eigenvalues from the negative ones and write
W k A ( T 4 A ( W W+ W.
Since W. 2 O, W > O, we can the Choleski decomposition write W+ O and
W- θ θ with + A_ Raxd Equation (10) follows then by direct computation: namely
(F) 2
1 Σ
f Sif
Σ Biy (f,, Wfz)
1
Σ
(fi, (Ω W)f. 2
Σ (fi, WF)
Σ
1 Σ J.j f, f Bzj (0,,, +
2
1
(
1 / Σ (f,, (Q W)()
1 Σ 0+
1 Σ - (
4
1
where we have used that Σig ā.
f 1+
B.2 ADDITIONAL. DETAILS ON LFD AND HED CHARACTERIZATIONS
In this subsection we provide further details and justifications for Definition 3.1 We first prove the following simple properties.
Lemma B.1. Assume we have a (continuous) process L + F(L) E na for t 2 O. The following equivalent characterizations hold:
(i) Dir (F())) + for (ii) EDr (F(L)IIF(III) + Pa/2 for + oo if and only if for any sequence - oo there exist subsequence ljx xo and unit limit F depending on the subsequence, such that Ar = Paff for r <d d
if and only if (t" + f.r ISrd.
Proof. Given F(t) E Raxd we can vectorize it and decompose in the orthonormal basis (e, @ of S r d, 0 kn ---1,, 1}, with {er canonical basis in Rd, and write
vee (())) Σ (t)e @ of cr,(() = (vec ((t)) @ 6
r,e
We then use Equation
25
to compute the Dirichlet energy
d n-
n-1
n-1
Dir ((
(1)x
- ΣΣ
(, Z gap A ΣΣ (t),),
r=1
where we have used the convention above that the eigenvector df is in the kernel of A. Therefore
n-1
(Fr)
(,((t) +0, t +c,
which occurs if and only if
(Ia Ø A)vec(F(t)) = ΣΣ R-1
Σ
(.) c, 8 Φ. +0.
17

---

Under review as a conference paper at ICLR 2023
(ii) The argument here is similar. Indeed write Q(t) [F(]] with Q(t) unit-norm signal. Namely, we can vectorize and write
vec Q())) = E E 0 Σ (1)) =1
Then & (Q()) - Pa/2 if and only if
( 7 DA, +
which holds if and only
Σ pr (f) +1
(1) -, E: <P.
(28)
given the unit norm constraint. This is equivalent to the Rayleigh quotient of Ig @ A converging to its maximal value Pa- When this occurs, for any sequence 1j extract converging subsequence that due to Equation (8) will converge to a unit eigenvector Q of La @ A satisfying (Ia @ A)Q contradiction that there exists sequence t, - oc such that 3 Ei (F(t, F/|[
have that (e,(()))) S
meaning that
Pa Q Conversely assume for Pa/2 c, for some > Then Equation (28) fails to be satisfied along the sequence, meaning that no subsequence converges a unit contradiction to our assumption.
eigenvector Foo of I A with associated eigenvalue Ps which is
Before we address the formulation of the system (t)) = AF in Raxd, with some initial condition F(0). We can vectorize the equation and solve véc (F())) [Ia @ (Fv))FF(t) meaning that
-omminatd dynamics, we solve explicitly
d カー
vec ΣΣ
e(1 (0)e, @ of
c,r(0) = (vec(F(0)), @ of
r=1 t=
Consider any initial condition F()) such that
Σ lo.ol Σ (vee @ o >0 0,
which satisfied for each vee vec(F(()) E add \ u-, where 4- is the orthogonal complement of
Rd & sp)n( Since LL- is lower-dimensional subspace, its complement dense. Accordingly for a.c. F()), we find that the solution satisfies
(Ivec(F(t)) 2²
Σ 2. + Ole (ee-)aa(a)))) - =c²²
+ O(e -
with P.l ( the projection onto Rd 2 ker(A), We that the norm of the solution increases exponentially, however the dominant term is given by the projection onto the lowest frequency signal and in fact
vee (Oc(F(M) (()))
PL A))(A
Veo veo((F())
+ (ap(e))p()))))) -
Keer ( ( )vec (F(O))
+ vec Foo):
vee F F(OII
(11
+ ((==s()p())))
such that (I. A Vec (Fo)) = 0 which means Af 0. for each column 1 1 d. Equivalently, one can compute Dir ((t)) and conclude that the latter quantity converges to zero as t + oo by the very same argument.
In fact, this motivates further the nomenclature LED and HFD. Without loss of generality we focus now the high-freqeency case. Assume that have HPD dynamics t. + F(I),
18

---

Under review as a conference paper at ICLR 2023
Le. ( /2, then vectorize the solution and write vec(F(tt) [[F(t]]vec(((() Q(t)), for some time-dependent unit vector rec(((t)) E Rnd:
vec (Q(()Q())))))) = Σ Ane(t)e, @ of
Σ (t) =
By Lemma B.1 and more explicitly Equation (28) we derive that the coefficients {ar,aa } associated with the eigevenctors e, & are dominant in the evolution hence justiffing the name high-frequency dominated dynamics.
The next result provides theoretical justification for the characterization of low (high) frequency dominated dynamics in Definition 3.1
Lemma B.2. Consider a dynamical system F(t) GNN (F((), with initial condition F(O).
(i) GNN, is LFD if and only if (Ia & A) vec(
[[F + 0 if and only if for each sequence t + do there exist subsequence tis - oc and Foo (depending on the subsequence) s.s. F(
[1 + Fo satisfying Af = 0, for each 1 r
(ii GNN is HFD if and only for each sequence and F (depending on the subsequence) s.t. [[ for each S S d.
+ there exist subsequence 1jx - F(tg)) + F satisfying Af PL fac.
Proof. (i) Since (1) - for each 1 r d if and only if (I, @ A)vec(F(1)) 0, we conclude that the dynamics is LFD if and only if (I& @ A) vec(F(()) - 0 due to (i) in Lemma B.1 Consider a sequence t, - co. Since (Fc()
[[F(t)
bounded sequence we can extract a converging subsequence t vec(F(t,, ()/(E (t,)) + Feo(F.o). If the dynamics is LFD, then vec(F(t;g)
(Ia 6 [[F III - O and hence we conclude that vee (oo) E ker(Ta @ A). Conversely, assume
that for any sequence I + oo there exists a subsequence 'js and Foo such that F(t;,)
satisfying Af. O, for each 1 S S d. If for contradiction had and t; [[P() II + Fx
such that eDir (F(L)//IFF()) 2 for large enough then by (i) in Lemma B.1 there exist S S d, > 0 and a subsequence tje satisfying
vec((((1))) er Φ 8(8) 00, [[((,))]
meaning that there is no subsequence of (4,, s.t. (Ia 60 Vee (1)1I) III + 0, providing contradiction.
(ii) This is equivalent to (ii) in Lemma B.I
Remark. We note that Lemma B.2 an L.F dynamics does not necessarily mean that the normalized solution converges to the kernel of @ A i.e. one in general has always to pass to subsequences. Indeed, we can consider the simple example + vec(F (t)) == cos(t) c, Φ for some S which satisfies (t) 0 for each but it is not convergent function due to its oscillatory nature. Same argument applies to HFD.
We will now show that Equation
36 can lead to a HFD dynamics. To this end, we assume that Ω W 0 so that Equation 36 becomes F(t) At) According to Equation (1)) the negative eigenvalues of W lead to repulsion. We show that the latter can induce HFD dynamics as
per Definition 3.1 We let Per be the orthogonal projection into the eigenspace of W & A associated with the eigenvalue p_ := [ I(ps 1). We recall that AW are the most positive and most negative cigenvalues of W respectively. We define EHFD by:
HFD min(p- AW I A), II W) 1)).
19

---

Under review as a conference paper at ICLR 2023
Theorem B.3. 1f NT then F(t) (t) is HFD for a.e. F(O):
E Ed (F((((()) =e 20
PA Pw P F(O)1 + Ole 2ICHED)
t20,
2
and F(L) [()Il converges Fo Raxd such that AP Pat for 1 S d. If instead P- < A. then the dynamics is LFD and ((t) converges to Fo E RRndd such that
A Ar for 1ss S d, exponentially fast.
Once we compute the spectrum of W & A via Equation 23 we can write the solution as recall that A In A so we can rephrase the eigenvalues A terms of the cigenvalues of A:
Proof of Theorem B.3
vec F())= = Σ (L (0)Φ OF
r.k
with W w & for S d, where (Φ,") )r is orthonormal basis of eigenvectors in Rd We can then calculate the Dirichlet energy along the solution as
EDrr (F(t)) 2 (vec (rt))F((t), (Ia @ A) vec F (t)))
Σ r,
e 22 11--
(0))
2
We now consider two cases:
If A > then X 12 ) S AW
If A. < O, then AW(1 (1- A²) S WA (Ps 1) = p-, with eigenvectors Φ. @ & for each W AW >" without loss of generality
assume that Pa is simple eigenvalue for A. In particular, if A, and A, (1 AA) p-, then
AL X max (A 1),
21(p 2
-1}}-
where W is the second most negative eigenvalue of and م is the second largest eigenvalue of A. In particular, we can write
n- PA gap(pat, - A), - aW pa W). (29)
From (i) and (ii) we derive that if A,W (1 - 1, # p-, then
x 2 p- < min(p- W P- W 11 1), 1A 11}}
min(p_ AW W gap(paI A), 8ap(1AW W) - 1)) -HHFD,
(30)
where we have used Equation
29
Accordingly, pf A then
(Di((F(()) = 26p-
PA
TA (0) +
Σ
20
2
r:A.
w
r.l:AW(1 (1-14)
21p-
PA
|PW F(0] + Ole (tteer)))))
By the same argument and for a.c. F(0) since Per (ec(F(()) = only vec(F(O)) belongs to a lower dimensional subspace of Rud;
factor out the dominant term and derive the following limit for t
vec(F (F(I) Pfr vec(F(0)) O(E cHP PP Ph )ve(((()))
pe vec (F(O))
vec (())
(1IPh v0)(F + ((e-2enrot)
w ved [(O)]]
where the latter is a unit vector voc F 2o) satisfying (I. & A)vec(F = PAvoc(F22)) which completes the proof.
For the opposite case the proof can be adapted without efforts as explicitly derived in the proof of Theorem 4.1
20

---

Under review as a conference paper at ICLR 2023
B.3 COMPARISON WITH CONTINUOUS GNNS: DETAILS AND PROOFS
Comparison with some continuous GNN models In contrast with Theorem 3.2 we show that three main linearized continuous GNN models are either smoothing more generally LPD. The linearized
PDE-GCN model Eliasof al. (021) corresponds choosing W and Ω W K(t) K(t))
in Equation 36). for some time- dene-dependent family t + K(t) E Rdxd;
F GCDE-GCNO()) = (1))(())(())))
The CGNN model Xhonneux al. (2020) be derived from Equation 36 by setting 52 I I- Ω -
=
F FGGNN(t) (t) F (t) F(t) + F(O).
Finally, in linearized BRAND(Chambert et al. (2021a) a row-stochastic matrix A(F(O)) is learned from the encoding via an attention mechanism and have
GTRAND()) nww (tt) (I ((F
We note that if A is not symmetric, then GRAND not a gradient flow.
Theorem B.4. PDE- GCND, CGNN and GRAND satisfy the following:
(i) PDE GCND is a smoothing model: & (FFDE-GON,(1) S
(ii) For a.e. F(0) it holds: CGNN never HFD and if we remove the source term, then Dir (FCGNN(t)||F CGNNNN)))) S -Rap( A)t
(iii) G is connected, FGRAND (t) - pl as oo, with fz" mean(f" (0)), S
By (ii) the source-free CGNN-evolution is LFD independent of Ω.. Moreover, by (iii), over-smoothing occurs for GRAND. On the other hand, Theorem 3.2 shows that the negative eigenvalues of W can make the source-free gradient flow in Equation (36) HFD. Experiments in Section 6 confirm that the gradient flow model outperforms CGNN and GRAND on heterophilic graphs
We prove the following result which covers Theorem 3.3
Proof of Theurem B.4 We structure the proof by following the numeration in the statement.
(i) From direct computation we find
di (F()) 2 di ((vec(F()) (La 8 A Avee
dt =
-vvvc((((t)) (K (t) @ A A vecc)(F(t))) I 0,
since K T (t)((()) & A 2 0. Note that we have used that (A @ BC D) AC & BD.
(ii) We consider the dynamical system
FCANN(() - AF() + (())+F(()))
We can write vec vec(F(()) E,, (0 x leading to the system
G.e(t) (1 ) cne(()), 0 D
We can solve explicitly the system as
caaf-> 1+
1
C.e(t) = c,)(
if if A
X
A
1
A
c,e(t) = c0, +t), otherwise.
We see now that for a.e. F(O) the projection (Ia 2 A
( T)vcc(F(tt) is never the dominant Pa. for any other non-maximal
Pa
term. In fact, if there exists s.t. X 2 Pa, then X A > X graph Laplacian eigenvalue. follows that there no Ω s.t. the normalized solution maximizes the Rayleigh quotient of I A, proving that CGNN is never HFD.
21

---

Under review as a conference paper at ICLR 2023
If we have source, then the CGNN equation becomes
F(t) - FF + F(t) = vec(F (t)) = ( @ (-A)))ee(((tt),
using the Kronecker sum notation in Equation 24 follows that we can write the vectorized solution in the basis {Øf & of Jr.e as
vcc(F(tt) = ci
Σ
G,0(0)² de + Ole sup(Afia-11y
Σ
(,,(0) @ &
ma
rap
Σ (0) 2 oه ) T.E>D
+ t
O(e ((E-AA))))))
meaning that the dominant term is given by the lowest frequency component and in fact, we normalize we find Dir (F(())))) S g-)a
(iii) Finally we consider the dynamical system induced by linear GRAND FGRAND(() - A F(t = -1 - (O))
Since we have no channel- mixing, without loss of generality we can assume that d = 1 - one can then extend the argument to any entry. We can use the Jordan form of A write the solution of the GRAND dynamical system as
f(t) = Pdiag(e"
dnt ) p-1r(0),
for some invertible matrix P of eigenvectors, with
m (me 1y
dit = (1-22
where my are the eigenvalue multiplicities. Since by assumption G is connected and augmented with self-loops, the row-stochastic attention matrix A computed in Chamberlain et al.) (2021a) with softmax activation is regular, meaning that there exists m N such that (Ar)) > for each entry (i.j). Accordingly, we can apply Perron Theorem derive that any eigenvalue of A has real part smaller than one except the eigenvalue Ag with multiplicity one, associated with the Perron eigenvector 1n. Accordingly, we find that each block ed'ut decays to zero as / -> oc with the exception of the one Jot associated with the Perron eigenvector. In particular, the projection of fo over the Perron eigenvector is just plln with y the average of the feature initial condition. This completes the proof.
B.4 PROPAGATING WITH THE LAPLACIAN
In this subsection we briefly review the special case of Equation (36) where Ω W, and comment on why generally expect framework where the propagation is governed by the graph vector field A to be more flexible than one with AA If 52 W and we suppress the source term i.e. W 0, the gradient flow in Equation 36 becomes
Ftt
(t
(31)
We note that once vectorized, the solution to the dynamical system can be written as
d
vec(F(()) = -Σ
(0) & Φ
In particular, we immediately deduce the following counterpart To Theorem 3.2
N Corollary B.5. If (ispe((W)))) R_ + 0, then Equation 331 HFD for a.e. F(O).
Differently from Equation (36 the lowest frequency component always preserved independent of the spectrum of W. This means that the system cannot learn eigenvalues of W to either magnify or suppress the low-fiee frequenegccc projection. In contrast, this can be done if Ω = O, or equivalently one replaces - A with A providing further justification in terms of the interaction between graph spectrum and channel-mixinn spectrum for why graph-convolutional models use the normalized adjacency rather than the Laplacian for propagating messages Kipf & Welling (007))
22

---

Under review as a conference paper at ICLR 2023
B.5 REVISITING THE CONNECTION WITH THE MANIFOLD CASE
In Equation 19 a constant nontrivial metric h in H leads to the mixing of the feature channels, We adapt this idea by considering a symmetric positive semi-definite H = W W with W E Rdxd and using it to generalize Dir by suitably weighting the norm of the edge gradients as
W F
Σ q
Σ
Σ
(VF")
1
Σ (A,))
IIW(F))
(32)
4
Ji(K.))EE
We note the analogy with Equation (9)) where the sum the nodes replaces the integration over the domain and the j-th derivative at some point i is replaced by the gradient along the edge (i,5) E. We generally treat W as learnable weights and study the gradient flow of ew
F(t) V (F(t)) -F())) W.
(33)
We see that Equation 33 generalizes Equation ① Proposition B.6. Let P be the projection onto ker ( W). Equation (33 is smoothing since
Ei (P))) S 2 W]gap(A)
|[(())] ( E ((PW @ ( 120.
In fact F(t) + Foo s.t. 3 x E R²: for each E V we have (foo): = vd + (.r((())
Proof of Proposition B.6
We can vectorize the gradient flow system in Equation (3) and use the
spectral characterization of W W & A in Equation (3 to write the solution explicitly
vec (ot)) E AP (0)) (0)Ø, W. & oF
where {AW), )r = spec (T)) W) C Rzu with associated basis of orthonormal eigenvectors given by (Ør W
h.. Then
EDi (F(t)) / (voc(F(tt), (Ia & VF Σ
2M 2 (O)E
1 Σ 2,0
(eee 20 XZ A
r: =0,e
r:A. >0
Dir ((Pwr & I veec +
Σ
(0)e 2M xP) X
2
r: >0,00
< Dir ((PW) @ (L))ve(F((O)) +
PA
21gap( W W)gap( (A) F(O))
2
where we recall that phrr is the projection onto ker W T W) and that by convention the index { = 0 is associated with the lowest graph frequency A = 0 0 by assumption G is connected. This proves that the dynamics is in fact smoothing. By the very same argument we find that
vec (t)) (I. & Pker vec() (F + W & L,)))c(F(()). t 00,
with pler the orthogonal projection onto ker the other terms decay exponentially to zero. We
first focus the first quantity, which we can write as
(Ia s pker vec (F(O)) > (.) s a
which has matrix representation of T E Rund with
do Σ _
By Equation (26 we deduce that the i-th row of of So E Buxd is the d-dimensional vector / Vi. We now focus on the term
(P & Ln)vec(F(0))
Σ
(r)(0), & o
23

---

Under review as a conference paper at ICLR 2023
which has matrix representation
Cr ) In particular, the i-th row is given by (,((0( 1,0, = PW (0).
P 0.
This completes the proof of Proposition B.6
Proposition B.6 B.6jimplies that no weight matrix W Equation (33 can separate the limit embeddings Fo. of nodes with same degree and same input features. In particular, we have the following characterization:
Projections of the edge gradients (VF) ( Rd into the eigenvectors of W W with positive eigenvalues shrink along the GNN and converge to zero exponentially fast as integration time (depth) increases.
Projections of the edge gradients (VF) (0) E IR into the kernel of W W stay invariant.
If W has a trivial kernel, then nodes with same degrees converge to the same representation and over-smoothingg Differently from Nt & Maehara (2019); Oono & Suzuki (2020): Cai & Wang (2020), over-smoothing occurs independently of the spectral radius of the "channel-mixing if its cigenvalues are positive even for equations which lead residual GNNs when discretized Chen et al. (2018).
According to Proposition [.6] we do not expect Equation 33 to succeed on heterophilic graphs where smoothing processes are generally harmful this confirmed in Figure (see prod-curve). To
deal with heterophily, one needs negative eigenvalues to generate repulsive forces among adjacent features.
A more general energy. Since in general one needs to generate repulsive forces too to deal with heterophilic graphs, we extend the Dirichlet energy associated with H W W > 0 to an energy accounting for mutual- - possibly repulsive interactions in feature space Rd, We first rewrite the energy eDir in Equation 32) as
Dir (F) = 2 E(E, (C,, W WTWWT WT) - 2
I Σ (C,, WTW))) i
(34)
If we replace the occurrences of W W with arbitrary symmetric matrices Z W E Rdxd we obtain
Es(F) :
1 Σ (f. S2F)
Σ
aty(P, Wf,) En (F) + cair
W (F).
(35)
2
with associated gradient flow of the form (see Appendix B
Ftt) (F(((F(t) - -()) + AF(!)W.
(36)
If we include the source term, then we have fully recovered the general energy in Equation c and its associated gradient flow.
C PROOFS AND ADDITIONAL DETAILS OF SECTION
We first explicitly report here the expansion of the discrete gradient flow in Equation 1I] after m layers to further highlight how this is not equivalent to single linear layer with a message passing matrix A as for SGCN Wu et al. 2019). For simplicity suppress the source term.
F(t +) = F(t) + (+) AF()tW)
(e +T) (Ind + (-
22 I2 W s A)) vec(F(L))
m
vec(F(mT)) = Σ m
& W A) & vee ((())
7*
(37)
and how the message passing matrix A actually enters the expansion after layers with each power 0 < & S m. This not surprising, given that we are discretizing a linear dynamical system, meaning that we approximating exponential matrix.
24

---

if and Under review as a conference paper at ICLR 2023
C.I FROM ENERGY TO EVOLUTION EQUATIONS: EXACT EXPANSION OF THE GNN SOLUTIONS
We first address the proof of the main result.
Proof of Theorem
4.1 We consider linear dynamical system
F(1 + F(1) AF())W, with W symmetric. We vectorize the system and rewrite it as
vec(F(t++)) (Ind TW & A)vec(F(tt)
which in particular leads to We can then write explicitly the solution
veo(F(mmT) (Ind W & veo
vec(F(mT))
(1+TA ^²))
(1) c,,(0) s of.
We now verify that by assumption in Equation (3) the dominant term of the solution is the projection into the eigenspace associated with the eigenvalue p_ | Kis 1). The following argument
follows same structure the proof of Theorem B.3 with the extra condition given by step-size. First, we note that for any such that A. > O, we have
+ |1+
M TA (1 x²
since required p- XW in Equation
13 Conversely, if A, 0, then
1 7 W (1
X² S max + TP-1, |,1 +>>
Assume that r|A > 1, otherwise there is nothing to prove. Then |1 + TP- > T|A only if
7 (2 pA) < 2,
which is precisely the right inequality in Equation [3] We can then argue exactly as in the proof of Theorem B.3 to derive that for each index such that A < 0 and A * AW, then
+ A x
max + + (pa 1)I, + 1 1)13
with AW and A defined in Equation 29
We can then introduce
2 _ (38)
6HFD = max p- - - Igap(pal A), - (Pa - ( W), AW
and conclude that
(m() = Σ (1- (1)
r,
(0n-1 oW +0 ((1+répr))
Σ
(,00 (1)o,)
(1-AM(1-)A)
=pp+T--)" (c)m W + +0
which completes the proof of Equation
Conversely, if p- W then the projection onto the eigenspace spanned by Φw @ @ is dominating the dynamics with exponential growth (1 +1+ (1+0))" We can then adapt the very same argument above by factoring out the dominating term once we note that due to the choice of symmetric normalized Laplacian A, we have &f Vd,, which then yields Equation [15
We also address the proof of Corollary
4.2
25

---

Under review as a conference paper at ICLR 2023
simply compute the Rayleigh quotient of Ia @ A. We report the explicit details for the HFD case since the argument for LFD extends without relevant modifications. Using Equation (12). compute the Dirichlet energy along a solution of F(L + r) = F(() - (() W satisfying Equation 13 by
Proof of Corollary 4.2 Once we have the node- wise expansion
E (F(m))))
Σ (1+ T XA)) (0)A
= (1 (1+ Tp- 2m
PA
ripa (0) +0
( TOHED) ((L+THEED)
11+ TP-
)) 2m
E
2
r:AW W
(1-A)(1))))
(1+TT__
PA ||PW F(O)112 O
Tp-
where Po is the orthogonal projector onto the eigenspace associated with the eigenvalue p- [AW Ips 1). In partieular, since
vec (((()))) = (1 + Tp-) (Pw vec(F(()) + O(8"))
we find that the dynamics HFD with vec(F()))/|locc(FF (t))l converging to the unit projection of the initial projection by Pw provided that such projection is not zero, which is satisfied for a.e. initial condition FO)).
C.2 COMPARISON WITH EXISTING RESULTS: PROOFS
Proof of Theorem 4.3 then
we drop the residual connection and simply consider F(t+T) = FFF(t)
vec(F(mr)) = (r W & vec (F()))
Since G not bipartite, the Laplacian spectral radius satisfies As < Therefore, for each pair of indices (r, E) we have the following bound:
1,, (1 X2 max
W 1)W),
The eigenvalues AW and AW are attained along the eigenvectors Φ W & Φ and Φw & Φ respectively. Accordingly, the dominant terms of the evolution lie in the kernel of L @ A, meaning that for any F, with non-zero projection in ker @ A) which is satisfied by all initial conditions except those belonging to a lower dimensional subspace the dynamics is LFD. In fact, without loss of generality assume that (W > X. then
and the inequality becomes strict if 0, i.c A2 >
vec (mmF)) = a W pm
Σ
PAW AW
PAW W T
((((m))
I
Σ
(Φ ( @ @dR)T voc(F(())
m:W==X
with (p(m) + 0 as m + oc, which completes the proof.
Gradient flow as spectral GNNs. We finally discuss Equation (I] from the perspective of
spectral GNNs as in Balcilar et al. (000)). Let assume that w 02 O. If we let A W A W (p ) be the eigendecomposition of the graph Laplacian and {AW] } be the spectrum of W with associated orthonormal basis of eigenvectors given by ( W }, and we introduce z" 2t() V R defined by 2/ (t) = (t), @W), then we can rewrite the discretized gradient flow as
z"(++ + w (I+ (I ) z' (()
z²(1) TA Az ((),
15r d. (39)
Accordingly, for each projection into the t-th eigenvector of W, we have a spectral function in the graph frequency domain given by A + 1 X4) If A. we have low-pass filter
26

---

Under review as a conference paper at ICLR 2023
while if A < O we have high-pass filter. Moreover,
that along the cigenvectors of W,
A, < then the dynamics is equivalent to flipping the sign of the edge weights, which offers
direct comparison with methods proposed in [Bo et al. (22IT) Yan et al. (202I) where some 'attentive' mechanism is proposed to learn negative edge weights based on feature information.
The previous equation simply follows from
(t 7) (f((LT)), = (f(() + W(Af()); 6,r)
21) W Σ (1),
which concludes the derivation of Equation 33
D PROOFS AND ADDITIONAL DETAILS OF SECTION 5
Proof of Theorem 5.7 First we check that time is continuous, then &g in Equation (35 is decreasing. We use the Kronecker product formalism to rewrite the gradient (FF)(F) as vector in R"d explicitly, we get
VE = (D I. W & A)voc(F) + (W @ (())).)(()).
It follows then that
(E((F(()) (VFE)(P(())) vec(F(())
dt
= (VFE)(F(t)) o (-V)Eo(F(tt))
If we introduce the notation Z(t) VFF((F(()) then rewrite the derivative
(F -Z(() (((1)((((())) (Z(()) = Σ (Z))) <
dt
by assumption on a. The discrete case follows similarly. Let the same notation above can write Fl. +) = F()) + (Z(()), with Z(() - &F(((F()))
Ea (F +)) (vec(F +)) ( I, w Ø A)vec(F(( + (W L.)vec(((O)))
1
( & I, - w
T)) (W 2 ve)cve(((()))
2
+ (o)()))
( I, - w w A) (vec (1) ra (W @ (.))ec(((())
1 (S2 @ I, w @ AF + (W @ Le)vee(F(()))
(vec(F(t)),
2
T(vec(F(tt)) ( I, w A)o(((t)))
+ + / ( 8 - w & () + (W 8 I.))ec(F(0)))
+ + (0 1. w & A)o (Z()))).
By using that Ω & I, W & is symmetric, we find that
= (() - W & A) + (W & L.)ve((((0)))
+ ((() (10 W A) (F( ()) -
= - ((t) + (F(t F(t), ( (n I W @ ())
1
-
S C+ F(O)11
where again have used that Z a (Z)) 0. This completes the proof.
27

---

Under review as a conference paper at ICLR 2023
To further support the principle that the effects induced by W are similar even in this non-linear setting, we consider a simplified scenario.
Lemma D.I. If we choose Ω = W = diag(w) with w" I for 1 r r S d and w = 0 i.e.
- (() solves the dynamical system
F(t) a (- (-AF(t)ddaa())))
with xx) 2 O, then the standard graph Dirichlet energy satisfies
de (E (F(())
0.
dt
Proof. This again simply follows from directly computing the derivative:
DE (F(t))
Σ
d
Σ (i,))E
fi(t)
50())
dt
Vd
vd,
r=1
d Σ Σ P= sEV
AP ( = Σ Σ y=l iEV
((W"")((WW((A))) 2 0.
Important consequence:
The previous Lemma implies that even with non-linear activations, negative cigenvalues of the channel-mixing induce repulsion and indeed the solution becomes less smooth as measured by the classical Dirichlet Energy increasing along the solution. Generalising this result to more arbitrary choices is not immediate and we reserve this for future work.
E ADDITIONAL DETAILS ON EXPERIMENTS
E.1 GENERAL EXPERIMENTAL DETAILS
GRAFF is implemented in PyToch Paszke et al. (019). using PyTorch geometric Fey & Lenssen (2019) and torchdiffeg Chen al. (018)). Code and instructions to reproduce the experiments are available on GitHub. Hyperparameters were tuned using wandt (2020) and random grid search. Experiments were run on AWS p2.8xlarge machines, each with Tesla V100-SXM2 GPUs
Methodology. Throughout the experiments and ablations we rely the following parameterisations. We implement v/N,, dDDE as single linear layers or MLPs, and we set Ω to be diagonal. For the real-world experiments we consider diagonally-dominann (DD) and diagonal (D) choices for the structure of W that offer explicit control over its spectrum. In the (Da)- we consider a W E dd symmetric with zero diagonal and E Fd defined by w. qa Eg W Vagl Ta- and set W diag(w) w) + W Due to the Gershgorin Theorem the eigenvalues of W belong to w. E,s Wo the model "can' easily re-distribute mass in the spectrum of [w,- I aasl
W via Ta- This generalizes the decomposition of W in Chen et al. (2020) providing a justification in terms of its spectrum. For (D) we take W to be diagonal. To investigate the role of the spectrum of
W on synthetic graphs, we construct three additional variants: W = W + W W = + = w named sum, prod and neg-prod respectively where prod (neg-prod) variants have only non nog-negative non-positive) eigenvalues.
E.2 ABLATION STUDIES ON SYNTHETIC HOMOPHILY
Synthetic experiments and ablation studies.
investigate our claims in a controlled environment use the synthetic Cora dataset of (Zhu al. 2020 Appendix G). Graphs generated for target levels of homophily via preferential attachment see Appendix E.3 for details. Figure confirms the spectral analysis and offers better understanding in terms of performance and smoothness of the predictions. Each curve except GCN - represents one version of W as in 'methodology" and we
implement Equation 11 with W = 0, Ω 0.
28

---

Under review as a conference paper at ICLR 2023
Figure 3 (top) reports the test accuracy vs true label homophily. Neg-pre Neg-prod is better than prod on low- homophily and viceversa on high-homophil
gm
This
confirms Theorem 4.1 where we have shown that the gradient flow can lead to HFD dynamics that are generally desirable with low-homophily through the negative eigenvalues of W. Conversely, the prod configuration (where we have an attraction-only dy- namics) struggles
homophily scenarios even
though residual connection is present. Both prod and neg-prod are "extreme' choices and serve the pur-
pose of highlighting that by turning off one side of the
spectrum this could be the more damaging depend- ing on the underlying homophily. In general though 'neutral' variants like sum and (DD) are indeed more flexible and better performing. In fact, (DD) outper- forms GCN especially in low-homophily scenarios,
Thee lube homaphly
Figure 3: Experiments on synthetic datasets with controlled homophily.
confirming Theorem 4.. where we have shown that without a residual connection convolutional models are LFD - and hence more sensitive to underlying homophily irrespectively of the spectrum of W. This is further confirmed in Figure
In Figure (bottom) we compute the homophily of the prediction (cross) for given method and we compare with the homophily (circle) of the prediction read from the encoding (i.c. graph-agnostic). The homophily here is a proxy to assess whether the evolution is smoothing, the goal being explaining the smoothness of the prediction via the spectrum of W as per our theoretical analysis. For neg-prod the homophily after the evolution is lower than that of the encoding, supporting the analysis that negative eigenvalues of W enhance figh-frequencies. The opposite behaviour occurs the case of prod and explains that in the low hhmoophlhll regime prod is under--erformant due to the prediction being smoother than the true homophily. (DD) and sum variants adapt better to the true homophily. We note how the encoding compensates when the dynamics can only either attract repulse (i.e. the spectrum of W has a sign) by decreasing or increasing the initial homophily respectively.
E.3 ADDITIONAL DETAILS ON SYNTHETIC ABLATION STUDIES:
The synthetic Cora dataset is provided by (Zuu et al. 2020 Appendix G). They use modified preferential attachment process to generate graphs for target levels of homophily. Nodes, edges and features are sampled from Cora propor- tional to mix of class compatibility and node degree result- ing in a graph with the required homophily and appropriate feature/label distribution. To validate the provided data be- fore use we provide Table 2 summarising the properties of the synthetic Cora dataset. All rows/levels of homophily have the same number of nodes (1,490), edges (5,936), fea-
100
3ym
Zym
Kom seen
0.4
2.5
08
19
tures (1,433) and classes (5).
Experiments on synthetic Cora Cora-GCN GCN ablation
Figure 4:
As well as the ablation shown in Figure 3 we used this dataset to perform an ablation using GCN as the baseline. We asses the impact of each of the steps necessary to augment standard GCN model to GRAFF. This involves 5 steps; 1) add an encoder/decoder. 2) add a residual connection. 3) share the weights of W and 1 across time/layers. 4) symmetrize W and Q2. 5) remove the non-linen -ineaiiyy between layers. The results are shown in Figure 4 annd corroborate Theorem 4.] that adding a residual term is beneficial especially in low-homophily scenarios. We also note augmentations 3,4 and 5 are not "costly" in terms of performance.
E.4 ADDITIONAL DETAILS ON REAL-WORLD ABLATION STUDIES
For the real-world experiments in Table we performed 10 repetitions over the splits taken from Pei ctt al. (2020). For all datasets we used the largest connected component (LCC) apart from Citeseer where the 5th and 6th split are LCC and others require the full dataset. For Chameleon and Squirrel
29

---

Under review as a conference paper at ICLR 2023
homophily
degree min we_regree
density
edge Sige_homophh
node_homoph
0.00
84.33
1.67
3.98
0.0027
0.00
0.00
0.10
71.33
2.00
3.98
0.0027
0.10
0.10
0.20
73.33
1.67
3.98
0.0027
0.20
0.20
0.30
70.00
2.00
3.98
0.0027
0.29
0.30
0.40
77.67
2.00
3.98
0.0027
0.39
0.39
0.50
76.33
2.00
3.98
0.0027
0.49
0.49
0.60
76.00
1.67
3.98
0.0027
0.59
0.60
0.70
67.67
2.00
3,98
0.0027
0.70
0.70
0.80
58.00
1.67
3.98
0.0027
0.78
0.79
0.90
58.00
1.67
3.98
0.0027
0.89
0.89
1.00
51.00
2.00
3.98
0.0027
1.00
1.00
Table 2: Summary of properties of synthetic Cora dataset
we added self loops and made the edges undirected as a preprocessing step. All other datasets provided undirected but without self loops. Each split uses 48/32/20 of nodes for training, validation and test set respectively. Table summarises each the datasets.
dataset
nodes
edges
features
classes
max degree
min degree
av degree
density
edge homoph
node homoph
Texas
183
558
1,703
5
104
3.05
0.0167
0.06
0.06
Wisconsin
251
900
1,703
122
3.59
0.0143
0.18
0.16
Cornell
183
554
1,703
5
94
3.03
0.0165
0.3
0.3
Film
7,600
53,318
932
5
1,303
7.02
0.0009
0.22
0.22
Squirrel
5,201
401,907
2,089
1,904
2
77.27
0.0149
0.23
0.29
Chameleon
2,277
65,019
2,325
5
733
2
28.55
0.0125
0.26
0.33
Citeseer *
3,327
9,104
3,703
99
0
2.74
0.0008
0.74
0.71
Citeseer
2,120
7,358
3,703
99
3.47
0.0016
0.73
0.71
Pubmed
19,717
88,648
500
3
171
4.5
0.0002
0.8
0.79
Cora
2,485
10,138
1,433
168
4.08
0.0016
0.8
0.81
Table 3: Summary of properties of real-word datasets. All LCC except
We perform an ablation study to further corroborate the behaviour seen in Figure For heterophilic datasets used the splits from Pei et al. (2020). For homophilic datasets used the methodology in Shhhhr et al. 2018). each split randomly selects 1,500 nodes for the development set, from the development 20 nodes for each class taken the training set, the remainder the validation set. The remaining nodes outside of the development set are used as the test set. This gives lower percentage (3-6%) of training nodes. This approach was taken because less training information is needed in the homophilic setting and performance can become less sensitive to other factors, meaning less signal from the controlled variable. We tested the structures of W against the real-world datasets with known homophily, again neg-prod outperforms prod in the heterophilic setting and vice-versa due the sign of their spectra.
allocated as
dataset
neg prod
prod
sum
Chameleon
67.32
58.86
68.36
Squirrei
51.39
42.11
51.29
Cora
31.80
79.65
81.17
Citeseer 32.47
67.31
67.53
Table 4: Ablation with controlled spectrum of W on real-world datasets
To validate the complexity analysis in Section 6 we performed a runtime ablation for the models between standard GCN and GRAFF described in the GCN ablation Figure 4 The average inference runtime over 100 runs for split of Cora was recorded. We also include runtimes for the provided dense and sparse implementations GGCN Yan et al. (202T). Adding the encoder/decoder (step 1) speeds up the model due to dimensionality reduction. Subsequent steps also reduce complexity and offer speedup with GRAFF performing the fastest.
30

---

Under review as a conference paper at ICLR 2023
0.5
0.4
OGCNL
0.2
0.0
Figure 5: Runtime ablation for inference on Cora dataset
E.5 DETAILS ON HYPERPARAMETERS
Using wandb Biewald (0200) we performed a random grid search with uniform sampling of the continuous variables. We provide hyperparameters that achieved the best results from the random grid search in Table An implementation that uses these hyperparameters is available in the provided code with hyperparameters provided in gra ff. www.ppaaaaspp... .y. Input dropout and dropout are rates applied to the encoder/decoder respectively with drupout applied in the ODE block. Further hyperparameters decide the use of non-linearities, batch normalisation, parameter vector w and source term multiplier / which are specified in the code.
w_style
Ir
decay
dropout
input input.dropon dropout
hidden
time
step_size
chameleon
diag .oo
0.0050
0.0005
0.36
0.48
64
3.33
squirrel
dimg_dom
0.0065
0.0009
0.17
0.35
128
2.87
texas
dim .dom
0.0041
0.0354
0.33
0.39
64
0.6
0.5
wisconsin
diag
0.0029
0.0318
0.37
0.37
64
2.1
0.5
cornell
diag
0.0021
0.0184
0.30
0.44
64
2.0
film
diag
0.0026
0.0130
0.48
0.42
64
1.5
Cora
diag
0.0026
0.0413
0.34
0.53
64
3.0
0.25
Citeseer
diag
0.0001
0.0274
0.22
0.51
64
2.0
0.5
Pubmed
diag
0.0039
0.0003
0.42
0.41
64
2.6
Table 5: Selected hyperparameters for real-world datasets
F.6 EXPERIMENTS ON LARGER HETEROPHILIC GRAPHS
Since our frameworks GRAFF and GRAFF are fast and efficient, we test them on two larger
heterophilic datasets arXiv-year Hu 2021 and snap-patents (Leskovec & Krevl 2014 Leskoved
et al. 2005)
datasct
nodes
edges
features
classes
edge homoph
arxiv-year 169,343
1,166,243
128
.222
snap maap-patents 2,923,922
13,975,788
269
5
.073
Table 6: Statistics of the larger heterophilic datasets.
Both the datasets concern predicting 'publication time" in a citation network. We report baselines as in Lim al. (202 and the best performing number among GRAFF and GRAFF in Table 7 Before commenting on the results, we observe that LINK (Zheleva Getoor 2009) is a method that only acts the input adjacency and ignores completely the node features. The facts that LINK is such a strong baseline on the datasets and that the MLP instead very poor denote that the features these datasets actually carry very little valuable information. Indeed, most of the MPNNs struggle
31

---

Under review as a conference paper at ICLR 2023
arXiv-year
snap-patents
MLP
36.7 + 0.21
31.34 = 0.05
L Prop 1-hop
43,42 + 0.17
30.28 + 0,09
L Prop 2-hop
46.07 = 0.15
38.61 0.07
LINK
53.97 # 0.18
60.39 = 0.07
SGC 1-hop
32.83
30.31 0.06
SGC 2-hop
3.06
29.09 0.09
C&S 1-hop
4.16
35.55
C&S 2-hop
19.78=0.22
49.08 0.04
GCN
46.02 0.26
45.65 + 0.04
GAT
46.05 0.51
45.37 0.44 (M)
GCNJK
46.28 0.29
46.88 # 0.13
GATJK
4.72
44.78
APPNP
38.15 = 0.26
32.19 + 0.07
H2GCN
49.09 + 0.1
(M)
MixHop
51.81 0.17
$2.16
GPR-GNN
45.07 0.21
40.19 # 0.03
GCNII
47.21 # 0.28
37.88 # 0.69
LINKX
56 + 1.34
61.95 + 0.12
GRAFF
50.52 = 0.23
55.77 = 0.27
Table Performance on larger heterophilic datasets. (M) stands for out of memory.
significantly, with the partial exception of MixHop which accounts directly for 2-hop information at each layer resulting in worse complexity that GRAFF. The reason why MixHop more competitive than other MPNNS is because these datasets have some strong 'monophily' type of bias such that architectures that directly access the more homophilic 2-hop are at an advantage; indeed, it has already been noted that LINK is indeed meant work well in this setting (Altenburger & Ugander] 2018). We also note that the labelling in arXiv-year was introduced in Lim et al. (2021). We point out how GRAFF manages to stay competitive with other MPNNs and almost consistently outperforms them. This confirms the validity of our theoretical analysis and motivates investigating energy functionals that would allow to incorporate higher-order terms as in MixHop in a gradient flow framework. We also note that by losing the inductive bias of message- pessaag-ppssing LINKX is not an optimal framework for dealing with homophilic graphs as reported [im et al. (021): on the other hand, an advantage of the gradient flow graph convolutional equations is that they manage to adapt to the underlying homophily of graph in a way that is provable and justifiable. The latter point is also in stark contrast with LINKX which instead works 'black box' on both adjacency and features.
Remark. Both arriv-year and snap-patents graphs come as directed and such information is essential for the task as already noted in Lim et al. (2021). Strictly speaking, by using the directed (and hence non symmetric) adjacency matrix A in Equation 17) and Equation (18) we are no longer gradient flow due the lack of symmetry of the gradient of the quadratic energy Eo noted in Appendix B.I It is still interesting to observe how symmetrizing and sharing the weights have enabled our framework to beat both GCN and GCNII by a margin even though are longer gradient flow. We reserve a more thorough investigation of gradient flows for directed non symmetric relations future work.
E.7 FURTHER COMPARISONS WITH BASELINES
In this subsection we provide further comparison with recent baselines (Luan et al. 2021 Lingam et al., 2021 Maurya et al.) 2021) that try to specifically target heterophily. We report their best numbers without reproducing them.
few comments are in order. First, we have reported the best numbers over all (several) configurations of the reported baselines. We note that despite its simplicity, framework is very competitive
32

---

Under review as a conference paper at ICLR 2023
Wiscomiin
Film
Squicrel 0.22
Chameleon
Hom level
0.21
0.30
0.22
0.23
0.74
0.80
0.81
#Nocks
183
251
183
7,600
5,201
2,277
3.327
18,717
2,708
#Edges
295
466
280
26,752
198,493
31,421
4,676
44,327
5,278
#Classes
5
77.15
HLP
4:22
TT.48
PSGNN
82.30
5.55
58.43
1.22
87,03 4: 5.77
36.07
0.69
73.48
2:13
78.14
1.25
77.19 1.35
89.73
0.39
81.73
1.36
GRAFF
85.38 1. 4.53
88.83
3.29
$1.05 1: 6.10
ar. 11 1.08
58.72
1 0.81
71.08
1.75
77.30
L8S
90.04
0.41
88.01 1: 1.00
GRAFFSL
86.40
84
87.26
1 2.52
77.50 3.24
35.96 6 0.95
59.01
1.51
71.38
1.47
76,81
L.12
89.81 1 0.50
87.81
1.13
Table 8: Node-classificatinn results.
Cora
Citeseer
Chameleon
Squirrel
Hom level
0.81
0.74
0.23
0.22
#Nodes
2,708
3,327
2,277
5,201
#Edges 5,278
4,676
31,421
198,493
#Classes
6
7
5
Test accuracy
87.69 + 0.86
76.75 + 1.33
70.44 + 1.47
59.60 + 1.49
Dir
0.41
0.22
0.94
0.80
E
Dir
0.16
0.10
0.73
0.77
1.93
1.02
0.98
1.02
A
-0.40
-0.95
-5.05
-4.80
Table Spectral analysis.
the small heterophilic graphs and stronger on the larger dataset Film and on all the homophilic ones. On the other hand, HLP and FSGNN are much stronger baselines on Squirrel and Chameleon. This also mainly because both the architectures handle the graphs as directed without self-loops which helps the performance massively and we suspect effectively reduces the heterophily of the graph hence making tha task easier. We reserve the investigation of heterophily in the context of directed graphs for future work.
E.8 SPECTRAL PROPERTIES OF THE LEARNT CHANNEL-MIXING
In Table 9 we report a few spectral properties of the channel-mixing W learnt on 4 real datasets two homophilic and two heterophilic along with the value of the normalized Dirichlet energy (X) |2 used for characterizing LFD and HFD dynamics.
X g (e)iii)//ll
Some important comments are order. The quantities Dir and &Dir are the normalized Dirichlet energy of the encoded node features before and after the diffusion layers respectively; recall that these are values between 0 and we introduced in Section 3 to characterize whether given GNN is mostly smoothing not. AW | instead is the ratio between most positive and most negative eigenvalue. (i) we note how the initial encoder provides us with node-wise features that are increasingly less smooth depending to adapt to the underlying homophily as we see from the fact that the normalized Dirichlet energy decreases much more after the message-passsingg for Cora and Citeseer. (iii) The ratio X / which is partial indicator of the amount of attraction vs repulsion exerted by the diffusion - is also indicator of performances
the underlying graph heterophily - as we derive from E (ii) The model learns
read when comparing the of CORA with the other datasets. (iv) As we increase the heterophily, the most negative eigenvalue increases in absolute value since most likely more repulsion is needed. We argue that the model might benefit from inducing fast attraction along some directions and fast repulsion along others since on heterophilic graphs we have larger positive and negative eigenvalues. Connected this last point though, we note that at the end of the architecture we have decoder that can also potentially discard some node feature projections hence keeping only the smoother (or less smooth) components.
33


Reviews:
[{"review_id": "eRVg2QoXiPv", "reviewer": "Reviewer_DfeD", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors present an approach of viewing the (residual) GNNs / graph ODEs through the point of an energy and corresponding gradient flow. In other words, the residual GNNs are viewed as derived from taking the gradient of a parameterised energy function. One of the ‘take-home’ points emphasised throughout the paper is the aspect that one should parameterise the energy and not the evolution equation of the GNNs.\n", "summary_of_the_review": "Below are a few comments that could hopefully improve the quality of the paper. They are mainly related to the intuition, motivation, and validation of the approach.\n\nMotivation of the energy view:\n\n- The energy functional itself is arrived in a sort of reverse-engineered manner by starting from the evolution equation (5), without the non-linearity $\\sigma$. Equation (5) itself seems to be presented without explaining why such a specific form alone would be necessary, but that does not affect the analysis seriously. Once the energy functional is defined, we are shown how the mixing matrix and its spectra regularizes the frequency spectra asymptotically — that such an energy helps mediate both smoothness (through positive eigenvalues of the mixing matrix W) and non-smoothness or high-frequency behaviour (through negative spectra of W), which is often connected to heterophily. This is an interesting concept, but not developed fully.  From Section 3 onwards, we see that the rest of the analysis takes the form of a regular graph/neural ODE such as GRAND, and does not directly make use of the energy except in Section 5.\n\n- In the end of Section 3, the authors argue ‘We argue that energies rather than evolution equations should be the object to parameterise for deriving more principled GNNs that are easier to interpret and analyse. ‘ This argument continues throughout without really being made clear of its necessity or advantages. \n\n- I like the perspective of the parameterised inspired from ‘finding the right notion of smoothness for our task’ as mentioned on page 4. However, beyond this minor intuition, I am unable to see the importance or the unique bias the energy view of GNNs brings over the evolution equations.\n\n- Lot of the later arguments regarding LFD/HFD is done using Dirichlect energy and seems disconnected from the energy functional we are introduced to before. This characterisation is nonetheless interesting and of value. However, the connection between HFD and heterophily is not obvious, particularly given that we are not shown how these properties pop up in the actual benchmark datasets, specially in low-homophily ones.\n\nExperimental validation \n\n- The authors present to us the interesting case of a synthesized Cora dataset at different levels of homophily. We are shown comparison with cases where the mixing matrix W is strictly negative in spectrum (neg-prod), positive in spectrum (‘prod’), and a general W. We are shown that in low homophile cases, the ‘neg-prod’ seems to have better performance over ‘prod’, and the other way round for high-homophily. However, in all cases, the individual performances never exceeds having a general W which means even in extremely low homophily situations, the energy functional learns to retain a component which promotes smoothness across features of LFD nature. This is interesting, and not entirely clear why - given the argument that heterophily essentially would be better off to have ‘repulsion’ amongst the particles.\n\n- The authors should better explain why the performance still remains lower than an MLP. Further, it would be advisable to show the reader the relative spectral components (signs of eigenvalues) of the learned W for different homophily levels — this would help validate the hypothesis that indeed the attraction repulsion dynamic is in play (as discussed in section 3.2). Also given that Graph isomorphism convolution layer (GINConv) [R0] is known to be one of the most expressive graph convolution layers available, I suggest to make comparisons with that also as a baseline in the experiments.\n\n- Similarly, it would be interesting to see the spectral aspects of the mixing matrices learnt for all of the real-world  heterophily datasets reported in Table 1 - do the learnt W and their eigenvalues mostly agree with the sign that is expected from the proposed analysis through Section 3.2?\n\n\nOther minor comments:\n\n- Theorem 5.1, I believe should be seen as Lipschitz continuity and not monotonicity as the authors claim on page 7.\n\n- The performance of GRAFF and GRAFF_{NL} is often fairly close — does this mean that the non-linearity  $\\sigma$ has relatively no role to play? Alternatively, how does the spectrum of the learnt W matrices differ for both the cases for a given dataset?\n\n- Before Section 5, the authors statement ‘Convolutional GNN models can deal with heterophily if the channel mixing matrix has negative eigenvalues. ‘ — does not seem clear or substantiated. The authors should make this clear, or explain this through examples.\n\n- In the experiments on benchmark data, we observe that in many cases, GGCN or even GCNII (that are known to be special cases of equation (5)) outperform GRAFF even at low homophily levels. Can this be explained in some way? Perhaps through the associated learnt energy functionals — currently I don’t see how the energy functional view has helped us in understanding these observations.\n\n\n- Finally, since the entire framework is based on dynamical systems, it would be nice to see if the model is able to capture actual dynamic systems, as considered in [R1] for example.\n\n\n\n[R0] How Powerful are Graph Neural Networks? K. Xu, W. Hu, J. Leskovec, S. Jegelka. ICLR 2019\n\n[R1] Chengxi Zang and Fei Wang. 2020. Neural Dynamics on Complex Networks. In Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining (KDD '20). Association for Computing Machinery, New York, NY, USA, 892–902. https://doi.org/10.1145/3394486.3403132\n", "strength_and_weaknesses": "Strength: The topic is definitely timely and interesting, and the general approach relevant. The paper is generally well-written. There are multiple concepts of interest, as using Dirchlect energy to quantify HFD/LFD nature of node features, emphasising the importance of residual connections, the attractive-repulsive multi-particle perspective. \n\nWeaknesses: The motivation, particularly with respect to the core idea of energy is unclear. Out of oversquashing, over-smoothing, and heterophily, the authors tell us that their approach is suited to arrive at principled GNNs, guided by the perspective of energy functional/gradient flow. However, the relevance of the energy function/ gradient flow to this is not entirely clear - in fact most of their analysis, proofs, and lemmas, continue to be based on the evolution function, and on the Dirichlect energy, that is used to analyse high-frequency and low-frequency asymptotic behaviour of the latent feature evolution.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "R5uE2lx3yWY", "reviewer": "Reviewer_Sead", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This mostly theoretical work presents a framework of GNNs in the context of Dirichlet energy minimization/maximization. GNNs are framed as estimators for gradients of energy functionals that point to either of those two directions. It shows that for certain classes of graphs and GNNs, attractive or repulsive forces between node features overweigh the other, leading to dominance of low or high frequencies, respectively, after convergence.  Based on this framework, the authors are able to design a GNN architecture that works well for low-homophily and high-homophily, depending on parameterization.", "summary_of_the_review": "In summary, I think this paper should be accepted, as it provides a new perspective of GNNs as gradient estimators for energy functionals. The paper, while reading, gave me a few ideas and a mental framework to work with, when thinking about GNN solutions for specific problems.\n\nQuestions:\n- This work focusses on architectures favoring low or high frequency content. I was wondering if it wouldn't be best to converge to a range of different frequency content, e.g. in different feature groups, to extract most from the given graph. What do the authors think here?\n- As far as I understand, the architecture only allows for switching between low and high frequency convergence with sign flipping. Could you think of a way of providing more explicit control over the frequency in the given framework?", "strength_and_weaknesses": "Strengths:\n- The presented theoretical framework is novel and interesting to read. It provides a new perspective on GNNs.\n- The framework unifies many existing GNN architectures and is able to derive conclusions about their behaviour with respect to favoring high or low frequencies. These conclusions are partially supported by proof of concept experiments.\n- The paper exposes a guideline of how to design GNNs that favor high frequency encodings, which might be useful to certain practical tasks, as most existing GNNs favor low frequency.\n- The paper is very well written. I am not an expert in ODEs and only roughly familiar with the most recent GNN theory and I was able to follow the paper mostly.\n\nWeaknesses:\n- The experiments are on the weak side of the paper. While theoretically shown effects are supported by proof of concept experiments, real applications of the presented phenomena are lacking.\n- In general, it would be great if the authors make stronger connections to practice, e.g. by zooming in on a specific application/graph type that would benefit from the presented method in a special way. That would also give a good intuition about how the knowledge can be applied to other tasks as well.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "QnNvmrzfzh", "reviewer": "Reviewer_Hwvc", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This work suggests an energy based perspective on GNNs, and by an analysis of the application of GNN layers and the channel-mixing operator sheds light on the behavior of GNNs.\nSeveral experiments are conducted showing improvement over the considered methods in some cases.", "summary_of_the_review": "The paper continues the interesting direction of interpreting GNNs through their energy function. I think that it has a merit but the experiments can be extended and compared with more methods, and many of the details are missing.", "strength_and_weaknesses": "Strengths:\n\n+This work takes inspiration from recent works like GRAND, PDE-GCN and GraphCON and extends the analysis of GNNs and their interpretation as dynamical systems that can be explored through their energies.\n\n+The authors sufficiently cite and discuss relevant background.\n\nWeaknesses:\n\n-The authors discuss the possibility of 'negative edge weight' in page 7. However, it is not clear what is the relation between allowing negative edge weights and recent works that propose designated architectures for this kind of data, for instance \"Learning signed network embedding via graph attention\" (AAAI 2020), \"Signed graph attention networks\" (ICANN 2019).\n\n-The authors show in Theorem 5.1 that the considered energy decreases from layer to layer which as discussed in the paper it is not anticipated to work well for heterophilic datasets. Then, I do not understand the obtained accuracy on such datasets. What causes the method to behave well? did the authors use only a small number of layers?\n\n-The phenomenon of oversmoothing is widely spread in GNNs and many works (that are properly cited in the paper) propose various ways to tackle this issue. The authors also discuss this issue but no study was done in this important aspect. \n\n-In page 8, paragraph \"The model and the parameterisation\" the authors suggest two models for the learned matrix W, but it is not clear what is the real influence of each of them. According to the appendix, I can only know what model was used, but as a reader I cannot infer what is the impact of those models on the different datasets. I think that the authors should add an experiment that studies this aspect.\n\n-It is not reported how many layers are used for the different experiments. I think that in this aspect it is important also given the theorems that regard to the energy change between layers.\n\n-Regarding the results in table 1, I believe that the authors are missing comparison with recent works with better performance. For instance \"Is Heterophily A Real Nightmare For Graph Neural Networks To Do Node Classification?\", \"Simple Truncated SVD based Model for Node Classification on Heterophilic Graphs\" and \"Improving Graph Neural Networks with Simple Architecture Design\" obtain significantly better accuracy on some of the datasets (e.g., Chameleon and Squirrel)\n\n-I appreciate the number of experiments but I think that to allow a broader comparison with many existing methods the authors should also report the accuracy on the semi-supervised Cora, Citeseer and Pubmed datasets.\n\n-It is claimed that the code and hyperparameters are provided but I cannot find them.\n\n*A question to the authors: can this method be used for graph classification? \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "5cySL7Rw3bS", "reviewer": "Reviewer_zvZU", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "Inspired by gradient flows as (solution curves to) differential equations that minimize an energy functional, the authors view GNNs as a gradient flow equation of a parametric energy. The authors show that in graph convolutional models (GCN), the positive/negative eigenvalues of the channel mixing matrix correspond to attractive/repulsive forces between adjacent features. They also demonstrate how the channel-mixing can learn to steer the dynamics towards low or high frequencies, which allows to deal with heterophilic graphs.\nThe authors provide a use case to show experimentally how  the gradient flow framework is more efficient than GCN, and achieves competitive performance on graph datasets of varying homophily often outperforming recent baselines specifically designed to target heterophily.", "summary_of_the_review": "The paper presents an  interesting (albeit not completely new) perspective to view GNNs from the lens of gradient flows (of a energy functional).  While I appreciate this new lens to view GNNs, there are a lot of implicit assumptions, confusing notions, vague definitions and various formulas that are not well explained or justified. As such, it makes decipher the mathematical statements and proofs hard to follow and ascertain. There are a number of ad hoc constructions. It is unclear the proposed framework really provides any deeper insight into how GNNs work.", "strength_and_weaknesses": "Strengths\n   + Provides an interesting (albeit not completely new) perspective to view GNNs from the lens of gradient flows (of a energy functional).\n\n    + Detailed proofs of the mathematical claims made in the paper\n\nWeaknesses:\n   - There are implicitly assumptions that are not clearly stated and some of the terminologies seem to be loosely used (see the next section for more details); \n   -  There also appear to be gaps between the \"theoretical framework\" and how it is applied in practice (with discrete datasets). ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject