Understanding convolution on graphs via energies
Authors: Francesco Di Giovanni
James Rowbottom
Benjamin Paul Chamberlain
Thomas Markovich
Michael M. Bronstein
Venue: Accepted by TMLR ()

Abstract: Graph Neural Networks (GNNs) typically operate by message-passing, where the state of a node is updated based on the information received from its neighbours. Most message-passing models act as graph convolutions, where features are mixed by a shared, linear transformation before being propagated over the edges. On node-classification tasks, graph convolutions have been shown to suffer from two limitations: poor performance on heterophilic graphs, and over-smoothing. It is common belief that both phenomena occur because such models behave as low-pass filters, meaning that the Dirichlet energy of the features decreases along the layers incurring a smoothing effect that ultimately makes features no longer distinguishable. In this work, we rigorously prove that simple graph-convolutional models can actually enhance high frequencies and even lead to an asymptotic behaviour we refer to as over-sharpening, opposite to over-smoothing. We do so by showing that linear graph convolutions with symmetric weights minimize a multi-particle energy that generalizes the Dirichlet energy; in this setting, the weight matrices induce edge-wise attraction (repulsion) through their positive (negative) eigenvalues, thereby controlling whether the features are being smoothed or sharpened. We also extend the analysis to non-linear GNNs, and demonstrate that some existing time-continuous GNNs are instead always dominated by the low frequencies. Finally, we validate our theoretical findings through ablations and real-world experiments.

Published in Transactions on Machine Learning Research (MM/YYY)
Understanding convolution on graphs via energies
Francesco Di Giovanni' University of Cambridge
Jd405@cam..cc uk
James Rowbottom' University of Cambridge
jr908@cam.acc uk
Benjamin P. Chamberlain Charm Therapeuties
Thomas Markovich Cash App
Michael M. Bronstein University of Oxford
Reviewed on OpenReview: https: openrevie.. net/ forum? id= ud=uFFuFPP
Abstract
Graph Neural Networks (GNNs) typically operate by message-passngg, where the state of node is updated based on the information received from its neighbours. Most message- message-passing models act as graph convolutions, where features are mixed by a shared, lincar transformation before being propagated over the edges. On node-clasification tasks, graph convolutions have been shown suffer from two limitations: poor performance on heterophilic graphs, and over-smoothing. It is common belief that both phenomena occur because such models behave as low-pass filters, meaning that the Dirichlet energy of the features decreases along the layers incurring smoothing effect that ultimately makes features no longer distinguishable. In this work, we rigorously prove that simple graph-convolutiona. models can actually enhance high frequencies and even lead to an asymptotic behaviour we refer to as over-sharpening, opposite to over-smoot,ing. We do so by showing that linear graph convolutions with symmetric weights minimize a multi-particle energy that generalizes the Dirichlet energy; in this setting, the weight matrices induce edgo- attraction (repulsion) through their positive (negative) eigenvalues, thereby controlling whether the features being smoothed or sharpened. We also extend the analysis to non-linear GNNs, and demonstrate that some existing time-continuous GNNs are instead always dominated by the low frequencies. Finally, we validate our theoretical findings through ablations and real-world experiments.
1 Introduction
GraphNeural Networks (GNNs) represent a popular class of neural networks operating on graphs (Sperduti, 1993; Goller Kuchler, 1996; Gori al., 2005; Scarsclli et al., 2008; Bruna ct al., 2014; Defferrard al., 2016). Most GNNs follow the message-passing paradigm (Gilmer et al., 2017), where node embeddings are computed recursively after collecting information from the 1-hop neighbours. Typically, Message Passing Neural Networks (MPNNs) implement convolution on graphs, where messages are first acted upon by lincar transformation (referred to as channell-miring), and are then propagated over the edges by a (normalized) adjacency matrix (Kipf & Welling, 2017; Hamilton et al., 2017; Xu et al., 2019; Bronstein et al., 2021).
While issucs such as bounded expressive power (Xu et al., 2019; Morris et al., 2019) and over-squashing (Alon & Yahav, 2021: Topping al., 2022; Di Giovanni et al., 2023) pertain to general MPNNs, graphh-onvolutionall MPNNs have also been shown suffer from additional problems more peculiar to node-classificttin tasks:

---

Published in Transactions on Machine Learning Research (MM/YYYY)
performance heterophilic graphs ie. those where adjacent nodes often have different labels and over-smoolhing of features in the limit of many layers.
Several works have gathered evidence that graph convolutions seem to struggle on some heterophilic graphs (Pei et al., 2020; Zhu et al., 2020; Platonov et al., 2023). common reasoning is that these MPNNs tend to enhance the low-frequency components of the features and discard the high-frequency ones, resulting in a smoothing effect which is detrimental on tasks where adjacent nodes have different labels (Nt & Maehara, 2019). Accordingly, common fix for this problem is allowing agned-message passing to revert the smoothing effect (Bo et al., 2021; Yan et al., 2021). Nonetheless, the property that graph convolutions act low-pase filters is actually only proven in very specialized scenarios, which leads to an important question:
Q.1 Are simple graph convolutions actually capable of enhancing the high-frequency components of the features without explicitly modelling signed message- message-passing?
A positive answer to Q.1 might imply that, in the limit of many layers, simple graph convolutions can avoid over-smoothing. soothhho occurs when node features become indistinguishable in the limit of many layers (Nt & Machara, 2019; Oono & Suzuki, 2020). In fact, Cai & Wang (2020); Bodnar et al. (2022); Rusch et al. (2022) argued that over-smoothing is defined by the Dirichlet energy Dir (Zhou & Schölkopf, 2005) decaying to zero (exponentially) as the depth increases. However, these theoretical works focus on a classical instance of GCN (Kipf Welling, 2017) with assumptions on the non-lineas activation and the singular values of the weight matrices, thereby leaving the following questions open:
Q.2 Can over-smoothing be avoided by simple graph convolutions? Tf so, do graph convolutions admit asymptotic behaviours other than over-smoothing, in the limit of many layers?
Contributions and outline. Both the poor performance on het- erophilic graphs and the issue of over-smoothing, can be characterized in terms of a fixed energy functional (the Dirichlet energy EDir decreas- ing along the features Ftt) computed by graph-coneooutional MPNNs. Namely, if & EDir (F(()) decreases w.r.t. the layer 1, then the features become smoother and hence likely not useful for scparating nodes heterophilic graphs, ult ultimately incurring over-smoothing. But does actu- ally Dir (F(t) decrease for all graph-convolutiona models? In this work we show that for many graph convolutions, a more general, parametric energy Eg, which recovers DDr as a special case, decreases along the features as we increase the number of layers. Since this energy can also induce repulsion along the edges, we are able provide affirmative answers to both Q.1 and Q.2 above, thereby showing that graph con- volutions are not bound to act as low-pass filters and over-smooth in the limit. Namely, our contributions are:
Figure 1: Actual gradient flow dy- namics: repulsive forces separate heterophilic labels.
In Section 4 we prove that large class of linear graph convolutions gradient flows, meaning that the features are evolved in the direction of steepest descent of some energy &a, as long as the weight (channel-mixing) matrices are symmetric. We study & and show that it provides a physical interpretation for graph convolutions as multi-particle dynamics where the weight matrices induce attraction (repulsion) on the edges, via their positive (negative) eigenvalues.
In light of this, we provide an affirmative answer to Q.1, proving that the class of gradient-flow MPNNs can enhance the high-frequency components of the ures (see Theorem 4.3 for the time-continuous case and Theorem 5.1 for the discrete setting). Our analysis suggests that gradient- MPNNs may be able to deal with heterophilic graphs, differently from other time-continuous MPNNs which instead can only induce smoothing of the features (Theorem 5.4), provided that a residual connection is available (Theorem 5.3).
2

---

Published in Transactions on Machine Learning Research (MM/YYY)
In Section 5 we also answer Q.2 positively by proving that in the limit of many layers, linear gradient-flow MPNNs admit two possible asymptotic behaviours. If the positive eigenvalues of the weight matrices are larger than the negative ones, then we incur over-smoothing; conversely, the phenomenon of over-sharpening occurs, where the features become dominated by the projection onto the eigenvector of the Laplacian associated with the highest frequency (Theorem 5.1). This result proves that asymptotic behaviours other than over-smoothing are possible, even for graph convolutions, depending on the interactions between the spectra of the graph Laplacian and of the weight matrices. To determine whether the underlying MPNN is sooothhigg/shappening the features, we also introduce measure based on the normalized Dirichlet energy, of independent interest.
In Section 6 we extend our theoretical analysis to graph convolutions with non-linear activations, showing that despite not being gradient-flows. the same parametric energy &g (F))) introduced above, is monotonically docrcasing with respect to t, thereby preserving the interpretation of the woight matrices as edge-wise potentials inducing attraction (repulsion) via their eigenvalues (Theorem 6.1).
Finally, in Section 7 we validate our theoretical findings through ablations and real-world experiments.
Reproducibility. Source code can be found al: http://gghabb.ccmm//oooboootooia
2 Graph convolutions and the Dirichlet energy
In this Section, we first review the mussagg-passingg paradigm, with emphasis on the class of convolational-type MPNNs. We then introduce the Dirichlet eneryy on graphs, a known functional which can be used to study the performance of graph convolutions on hetcrophilic tasks and the over-smoothing phenomenon.
2.1 Preliminaries on graphs and the MPNN class
Notations and conventions. Consider graph G EV, with n := [V] nodes and E C CV V representing the edges. We assume that G is simple, undirected, and connected. Its adjacency matrix A is defined as agy = 1 if (i.j) E and zero otherwise. We let D diag(do. (do.dd....... do.....,n_.1) be the degree matrix. We typically use A to denote some (normalized) version of the adjacency matrix. We are interested in problems where the graph has node features {f; E Id ; E V} whose matrix representation is F E Rnxdd In particular, f; E IRd is the i-th row (transposed) of F, while f" E R" is ils r-th column. Later, we also rely on the vectorization of the feature matrix vec(F) E gnd simply obtained by stacking all colurns of F.
MPNNs and the 'convolutional flavour: In this work investigate the dynamics of MPNN in terms of its smoothing sharpening effects on the features, which properties particularly relevant in the context of heterophily and over-smoothing. Accordingly, we focus on node-level tasks, where cach node i E V has a label y/ E (I,...,K') K} that we need to predict. Message Passing Neural Networks (MPNNs) (Gilmer et al., 2017) compute node embeddings (t) at each node and layer t, using the following recursion:
(t + 1) = UP, (f,(), AGG, (L)()): : (1,5) E EB)),
where AGG, is invariant to permutations. The svor-smog phenomenon was studied for a class of MPNNs that we refer to gmaph convolutions (Bronstein et al., 2021), which essentially consist of two operations: applying shared linear transformation to the features ('channel mixing') " and propagating them along the edges ("diffusion"). Namely, we consider MPNNs whose layer update can be written in matricial form as:
F F(!) +
( (-F(t)&, + AF()),--((O))
where Zr W., and W, are learnable matrices in dxd performing channel mixing, while is pointwise nonlincarity; the message- paaaggg-ppsnnggggoo matrix A instcad, Icads to the aggregation of features from adjacent nodes. We first. note how this system of equations include common instances of the MPNN class with residual connections. If A = D 11/2 1/2 and Ω = W. = 0, then we recover GCN (Kipf & Welling, 2017). The case of A D A and W.- = 0 results in GraphSAGE (Hamilton et al., 2017), while by choosing Ω = o and W, and W, as convex combinations with the identity we recover GCNII (Chen et al., 2020). Finally, if A A, S -(1 + W, and W, 0, then (1) is a shallow variant of GIN (Xu et al., 2019).
3

---

Published in Transactions on Machine Learning Research (MM/YYN
2.2 The Dirichlet energy graphs
To assess whether graph convolutions induce smoothing or sharpening behaviour over the features, we monitor their Dirichlet energy. Consider a (feature) signal F - R4, where we associate a vector f; to each node E V. The graph Dirichlet energy F is defined by (Zhou & Schölkopf, 2005):
f
EDi
E (R,))E
I(VF),
(VF),
dd
where VF is the gradient of F and is defined edge-wise. We see that Dir measures the smoothness of F; meaning whether F has similar values across adjacent nodes.
The graph Laplacian. The (normalized) graph Laplacian is the operator - The Laplacian is symmetric and positive semidefinite, and its eigenvalues are referred to as (graph) frequencies and written in ascending order as λo X S... S An-1. It is known that À = 0 while An-1 S 2 (Chung & Graham, 1997). The low frequencies are associated with macroscopic (coarse) information in the graph, while the high frequencies generally describe microscopic (fine-grained) behaviour. In fact, if let {Φe} } be orthonormal frame of eigenvectors for A (i.e. A Xeσe), then features F with vector representation (ec E Rnd obtained by stacking all columns, as
rewrite the Dirichlet energy of the
E Dir (F) = -E
Xe f22= trace(F' AF) = (vec (F)) FI) & A)vec(F),
(2)
1=1 t=0
where & denotes the Kronecker product some properties are reviewed in Appendix A.2 and La is the identity matrix in Raxd, Feat are smooth if the Dirichlet energy is small, which occurs when for each channel r, we have large projections onto the eigenvectors Φe associated with lower frequencies.
3 Heterophily and over-smoothing: overview on related work
On a broad scale, our work related studying GNNs as filters (Defferrard al., 2016; Bresson & Laurent, 2017: Hammond et al., 2019; Balcilar et al., 2020; He et al., 2021) and adopts techniques similar to Cai & Wang (2020). Our approach is also inspired by Haber & Ruthotto (2018); Chen et al. (2018); Biloš et al. (2021), where layers of an architecture are regarded as Euler discretizations of time-contimuous dynamical system. This direction has been studied in the context of GNNs in several flavours (Xhonneux et al., 2020; Zang & Wang, 2020; Chamberlain et al., 2021a; Eliasof et al., 2021; Chamberlain et al., 2021b; Bodnar et al., 2022; Rusch et al., 2022); we share this perspective our work and expand the connection, by studying the energy associated with GNNs and think of features as 'particles' in some generalized position space.
Understanding heterophily via the Dirichlet energy. A nodo-clssification task graph G is helerophilic when adjacent nodes often have different labels. In this case the label signal has high Dirichlet energy, meaning that enhancing the high-frequency components of the features could be useful to separate classes - the prototypical example being the eigenvector of A with largest frequency An-1, that can separate the clusters of a bipartite graph. While MPNNs as in (1) have shown strong empirical performance on node-classification tasks with low heterophily, has been observed that simple convolutions on graphs could struggle in the case of high heterophily (Pei et al., 2020; Zhu et al., 2020). Such findings have spurred a plethora of methods that modify standard MPNNs to make them more effective heterophilic graphs, by by introducing some form of signed pesssae-passingg to
either ignoring rearranging the graph structure, 'reverse' the dynamics and magnify the high frequencies (Bo et al., 2021; Luan et al., 2021; Lim et al., 2021; Maurya et al., 2021; Zhang et al., 2021; Wang et al., 2022; Luan et al., 2022; Eliasof et al., 2023).
Understanding over-smoothing via the Dirichlet energy. Whether a given MPNN acts as a low-pass filter or not, also at the heart of the over-smoothing phenomenon, regarded as the tendency of node features to approach the same value- up degree rescaling as the number of layers increase, independent of any input information (Li et al., 2018: Oono & Suzuki, 2020). Cai & Wang (2020) formalized over-smoothing
4

---

Published in Transactions on Machine Learning Research (MM/YYYY)
the Dirichlet energy of the features E Eiir (F(t)) decaying to zero as the number of layers t diverge. A similar approach has also been studied in Zhou et al. (2021); Bodnar et al. (2022); Rusch et al. (2022) - for review see Rusch et al. (2023). Although over-smoothing is now regarded as a general plague for MPNNs, it has in fact only been proven to occur for special instances of graph convolutional equations that (i) have nonlinear activation ReLU, (ii) have small weight matrices (as measured by their singular values), and (iii) have no residual connections.
The general strategy. Both the problem of graph convolutions struggling on heterophilic tasks and of over-smoothing, arise when a fixed energy functional, in this case the Dirichlet energy, decreases along the features computed at cach layer. Namely, both these phenomena occur when, in general, d E (F(t)) S O, with F(L) the features at layer L. To address the questions Q.1 and Q.2 in Section 1, we sludy whether there exist general energy functionals promoting dynamics more expressive than simple smoothing of the features, that decrense along graph convolutions as in (1). Analyzing complicated dynamical systems through monotonieity properties of energy functionals is common in physics and geometry. Indeed, in this work we show that the same approach sheds light on the role of the channel-mixing weight matrices for graph-convolutions, and on the dynamics that MPNNs can generate beyond
de
4 Gradient-flow MPNNs: understanding convolution on graphs via parametric energies
In this Section we study when graph-convolutions (1) admit an energy functional &g such that t + (F())) is decreasing as L increases. This allows us to interpret the underlying MPNN-dynamics as forces that dissipate (or in fact, minimize) the energy Eg. By analysing Eg we can then derive qualitative properties of the MPNN, and determine if the features F(t) incur over-smoothing or not when t is large. Put differently, we investigate when MPNNs as in (1) are gradient flows, special class of dynamical systems we briefly review next.
What is gradient flow? Consider an N N-dimensional dynamical system governed by the differential equation F(t) (t)( that evolves some input F(0) for time f 2 0. We say that the cvolution equation is a gradient. flow if there exists E R-N - R such that F(F(()) = (V)(F()). In this case, since ((((t)) = -I (FF((t) the energy E decreases along the solution Ftt). Gradient flows are valuable since the existence of and the knowledge of its functional expression allow for a better understanding of the underlying dynamics. To motivate our analysis of graph convolutions, we first review examples of gradient flows graphs that do not involve learnable weights we discuss variational methods for image processing (You et al., 1996; Kimmel et al., 1997) in Appendix A.1 and Appendix B.5.
4.1 Gradient flows on graphs: the non-leareable case
A prototypical gradient flow: heat equation. Let F E Ruxd be the feature matrix. The most common form of difluusio process on G is the heat equation, which consists of the system (') = AP) where fr E F" is the r-th entry of the features for 1 S r S d. The heat equation is an example of gradient flow: if we stack the columns of into vec(F) E Rnn, we can rewrite the heat equation as
1 ( V (eo = - 2
(3)
where Dir Rnd + R is the (graph) Dirichlet energy defined in (2). The evolution of F(t) by the heat equation (3) decreases the Dirichlet energy &Dir (F(()) and is hence smoothing process; in the limit, eir (Di((F(t) + 0, which is attained by the projection of the initial state F(O) onto the null space of the Laplacian.
Label propagation. Given a graph G and labels (y:) on V1 C V, assume we want to predict the labels on V2 V. Zhou & Schölkopf (2005) introduced label propagation (LP) solve this task. First, the input labels are extended by setting (0) 0 for each i V \ V,. The labels are then updated recursively according to the following rule, that is derived as gradient flow of an energy ELP we want to be minimized:
1
(t) = - (Y()), &L. ELP(Y): =
: EDir
+ Y - Y(O)
(4)
5

---

Published in Transactions on Machine Learning Research (MM/YYYY)
The gradient flow ensures that the prediction is attained by the minimizer of E..P, which both enforces smoothness via Dir and penalizes deviations from the available labels (soft boundary conditions).
Motivations. Our goal amounts to extending the gradient flow formalism from the parameter-free case to deep learning setting which includes features. We investigate when graph convolutions as in (1) admit (parametric) energy functionals that decrease along the features computed at each layer.
4.2 Gradient flows on graphs: the learnable case
In the spirit of Haber & Ruthotto (2018); Chen et al. (2018), we regard the family of MPNNs introduced in (1) as the Euler discretization of a continuous dynamical system. More specifically, we can introduce step size T (0, 1] and rewrite the class of MPNNs
F(I+ F(() (-F()) (F()) F())W,),
which then becomes the discretization of the system of differential equations
F(t) = (-F()) + AF(!) F())..),
where denotes the time derivative. Understanding the mechanism underlying this dynamical system might shed light on the associaled, discrete family of MPNNs. By this parallelism, we study both the time-continuous and time-discrete (i.e. layers) settings, and start by focusing on the former for the rest the section.
4.2.1 Which energies are being minimized along graph coppolconvolttional models?
Assume that A in (1) is symmetric and let F(0) be the input features. We introduce a parametric function Eg : Rnd + IR, parameterised by x d weight matrices Ω and W of the form:
(F) E Stf.)
Ac WEJ) + ((F,)()))
(5)
goaic
Srouten
Note that we can recover the energies in Section 4.1. If = W = = and 60 = O, then Eg = Eir as per (2), while if pp is an Le-penalty, then &g = ELP as per (4). We can also recover harmonic energies on graphs
(see Appendix B.5). Importantly, we choose 40 (F, F(0)) = 2,, (f, WF,(0),, for W e Rddd then we can rewrite (see Appendix B.1)
(F) (vec(F), ( L, W & FF) + 2(W L,)vec(F(())).
(6)
The gradient flow of Ea then be simply derived as (once divide the gradient by 2):
1 F(() - 1100=---reee V FF()) Ω S2 (++2) + AF())
(W+
W WW+W W
FOO)..
(7)
Since Ω, W appear in (7) in symmetrized way, without loss of generality assume Ω and W to be symmetric d x d channel mixing matrices. Therefore, (7) simplifies as
F(t) - F(t AF(t)W- F(O)W.
(8)
Proposition 4.1. Assume that G has a non-trivial edge. The linear, lime-continuous MPNNs of the form
((t) = - (P(t))-F(0)))
are gradient flows of the energy (6) if and only if the weight matrices and W are symmetric.
Therefore, large class of lincar MPNNs evolve the features the direction of steepest descent of an energy, provided that the weight matrices are symmetric. Note that despite reducing the degrees of freedom, the symmetry of the weight matrices does not diminish their power (Hu et al., 2019).

---

Published in Transactions on Machine Learning Research (MM/YYY)
4.2.2 Attraction and repulsion: physics-inspiree framework
Thanks to the existence of an energy Eg, we can provide a simple explanation for the dynamics induced by gradient-flow MPNNs as pairwise forces acting among adjacent features and generating attraction and repulsion depending on the eigenvalues of the weight matrices.
Why gradient flows? A multi-particle point of view. Below, we think of node features as particles in Rd with energy Eg. In (5), the first term E is independent of the pairwise interactions and hence represents an 'external' energy in the feature space. The second term & instead accounts for pairwise interactions along edges via the symmetric matrix W and hence represents an energy associated with the graph st succuue For simplicity, we set the source term 80 to zero and write W = = O θ + O_O_. by decomposing it into components with positive and negative eigenvalues. We can then rewrite Ey in (5) as
EO(F) = (f., ()))
1 Σ (V+()) + 2 2
2 | (VF)
Σ i
(9)
graph-independent
adtraction
repulsion
which we have derived in Appendix B. Consider now the gradient flow (8) where features F(t) are evolved in the direction of steepest descent of Eg. Recall that the edge gradient (VF(t)), is defined edge-wise and measures the difference between features ((()) and (1) along (i,j) E. We note that:
(i) The channel-mixing W encodes allractive edge-wise interactions via its positive eigenvalues since the gradient terms |++ (VF(t))) are being minimized along (8), resulting in a smoothing effect where the edge-gradient is shrinking along the eigenvectors of W associated with positive eigenvalues;
(ii) The channel-mixing W encodes repulsive edge-wise interactions via its negative eigenvalues since the gradient terms O (VF((t)) are being minimized along (8), resulting in a sharpening effect where the edge-gradient is expanding along the cigenvectors of W associated with negative cigenvalues.
Next, we formalize the smoothing vs sharpening effects discussed in (i) and (ii) in a more formal way.
4.3 Low vs high frequency dominant dynamics: a new measure
Attractive forces reduce the edge gradients and are associated with smoothing effects which magnify low frequencies, while repulsive forces increase the edge gradients and hence afford a sharpening action enhancing the high frequencies. To determine which frequency is dominating the dynamics, we propose to monitor the normalized Dirichlet energy: (F)))(tt)) [())) This the Rayleigh quotient of I & A and so it satisfics 0 < Dir (F)/1[FI2 S An-1 (sce Appendix A.2). If this Rayleigh quotient is approaching its minimum, then the lowest. frequeney component is dominating, whereas if it is approaching its maximum, then the dynamics is dominated by the highest frequencies. This allows us to introduce the following characterization, of independent interest, to measure the reesooncyy-rsspons of spatial MPNNs:
Definition 4.2. Given a time-continuous (discrete) MPNN computing features F(t) at cach time (layer) t, we say that the MPNN is Lou-Prequenccy-DDominnw (LFD) if gir ( (F()) - for I + 0o. Conversely, (HFD) if Dir (F(t)/ + An-1 as + xc.
we say that the MPNN is
Note that the LFD characterization differs from the notion of over-smoothing introduced Cai & Wang (2020); Rusch et al. (2022), since it also accounts for the norm of the features. In fact, in Appendix B.2 we derive how our formulation also captures those cases where the Dirichlet energy is not converging to zero, yet the lowest frequency component is growing the fastest as time increases. Our notion of HFDD-yyaamcss is also novel since previous works typically focused on showing when graph convolutions behave as low-pass filters.
Similarly to previous works, we now focus on a class of graph-coneolutions inspired by GCN and show that even if we remove the terms Q2 and W from (I), caaph-convolutions can enhance the high frequencies and ultimately lead a HFD dynamics opposite to over-smoothing.
Assumption. We let A = D-122 AD-ADD-11 1// and consider the simplified gradient flows with Ω W 0.
7

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Theorem 4.3. Given a continuous MPNN of the form F(t) = AF()) let pm M S ... Ad-1 be the eigenvalues of W. Tf ( - 1) > Pa-1, then for almost every F(O), the MPNN is HFD. Conversely, if 1401(A2-1-)) - 1) < pd-1, then for almost every input F(O), the MPNN is LFD.
We provide convergence rales in Theorem B.3 in the Appendix. In accordance with (9), Theorem 4.3 shows that the channel-mixing matrix W can generate both attractive and repulsive forces along edges, leading to cither an LFD or a HFD dynamics. In fact, the MPNN is HFD when W has a negative cigenvalue Ao sufficiently larger than the most positive one pd-1 and viceversa for the LFD case. Therefore, Theorem 4.3 provides affirmative answers to Q.1, Q.2 in Section in the time-continuous case, by showing that a simple class of gradient flow convolutions on graphs can learn to enhance the high frequencies and lead to a HFD dynamics where the highest frequency components dominate in the limit of many layers.
We note that Eg can generally be negative and unbounded. With slight abuse of nomenclature, we call them 'energy' since the associated MPNNs follow the direction of steepest descent of such functional.
Remark.
5 The interactions between the graph and channel-mixing spectra
In this Section we consider discretized gradient flows, so that we can extend the analysis to the case of layers of an architecture. We consider MPNNs that we derive by taking the Euler discretization of the gradient flow equations in Theorem 4.3; given a step size T and W E Rddd symmetric, we have
F(t+T) F(t) + AFF()W, F(O) VPN (Fo),
(10)
where encoder VENN RR + Raxd processes input features Fo and the prediction (DD((((T) is produced by a decoder VDDE : Roxd_ + B Here, K is the number of label classes, T = mr is the integration time, and m is the number of layers. We note that (i) typically VEN, DE are MLPs, making the entire framework in (10) non-linear; (ii) since we have a residual connection, this is not equivalent to collapsing the dynamics into a single layer with aggregation matrix A Am as done in Wu et al. (2019) see (31) in the Appendix.
5.1 Discrete gradient flows and spectral analysis
The gradient flow in (10) is a linear, residual GCN with nonlincar encoder and decoder operations. Once we vectorize the features F(1) + ve)))(()) E Ren we can rewrite the update as (e = ve)(F(1)) + (Wa & A) vt)(F(t)) (sce Appendix A.2 for details). In particular, if we pick bases {4/) C Rd and (Φ.) C R. of orthonormal eigenvectors for W and A respectively, we can write the features after m layers explicitly.
d-1a-1
vec(F(mTT)
ΣΣ (1 (1+T4(0 -x)) (O)a0) @ Φb.
(11)
where (0) := (ve)(F(O)), bb, & 6i) and {ur} are the eigenvalues of W. We see that the interaction of the spectra {4r} and {ee] is the 'driving' factor for the dynamics, with positive (negative) eigenvalues of W magnifying the frequencies d < (> respectively). In fact, note that the projection of the features onto the kernel of W stay invariont. In the following we let no S #1 $1 ... S pa-1 be in ascending order. Note that Φn-1 is the Laplacian cigenvector associated with largest frequency An-1- We formulate the result below in terms of the following:
Pd-1
|Pm|l
2
(12)
An-1
7(2 112-xum -2--)
Note that if (12) holds, then 0①0 < 0 since An-1 1 < 1 whenever G is not bipartite. 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 approximate the gradient flow in (8). We restrict to the case where Po, Pa-1 are simple eigenvalues, but extending the results to the degenerate case is straightforward.
constraint on the step-size since if r too large, then we no longer
8

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Theorem 5.1. Consider MPNN with update rule F(t + r) = + 7FF(t)W with W symmetric. Given m layers, if (12) holds, then there erists ō < 1 8.1. for all E V we have:
+
(((T))= = 1100l000-----1111
(co.n-1 an up + 0 ))
(13)
Conversely, if Pa-1 > (pol(An-1 1), then
di (ct-1o(0)y (ca 2|E| +0-1+0
f((T)) (1 (1+TTag =Ma
(14)
We first note that 5 is reported in (32) in Appendix C.1. II (13) holds, then repulsive forces (high frequencies) dominate since for all E V, we have f((T) ~ 1 20 up to lower order terms in the number of layers. Thus as we increase the depth, any feature (((r)) becomes dominated by a multiple of the E Rd that only changes based on the value of the Laplacian eigenvector Φn-1 node i. Conversely, if (14) holds, then (((T)) ~ vd:- - meaning that the features become dominated by multiple of a only depending the degree of which recovers the over-smoothing phenomenon (Oono & Suzuki, 2020).
Corollary 5.2. If (13) holds, then the MPNN is HFD for almost every F(O) and F to Foo s.t. Af = >n-1ff for each r. Conversely, if (14) holds, then the MPNN is LFD for almost every F(() and (((mr)/[[((mr) converges to Foo s.l. AP 0 for each r.
converyes
Over-smoothing and over-sharpening. Our analysis shows that (i) linear grapl-connolltional equations can induce a sharpening effect (13) (which answers Q.1 in Section affirmatively); (ii) graph convolutions can avoid smothhng through the negative eigenvalues of the channel-mixing by incurring the opposite behaviour of over-sharpening (i.c. HFD behaviour) in the limit of many layers, which provides affirmative answer to Q.2. As per Corollary 5.2, over-sharpening entails that the features converge, up to rescaling, to the eigenvector of A associated with the largest frequency. Similarly, Corollary 5.2 implies that LFD is a rigorous charactcrization of over-smoothing, given that features converge to a vector with zero Dirichlet energy in the limit of many layers. Accordingly, both over-smoothing and over-sharpening afford a loss of information, since the features converge to the Laplacian eigenvector with minimal or maximal Dirichlet. energy, respectively. When the number of layers is small though, the features F(mr) also depend on the lower-order terms of the asymptotic expansion in Theorem 5.1; whether the dynamics is LFD or HFD will then affect if the lower or higher frequencies have larger contribution to the prediction.
The role of the residual connection. The following result shows that the residual connection is crucial for the emergence of the over-sharpening (HFD) regime:
Theorem 5.3. G is not bipartite, and we remove the residual connection, iee. F(I + W, with W symmetric, then the dynamics is LFD for almost every F()), independent of the spectrum of W.
We see that the residual connection enables the channel-miaing to steer the evolution towards low or high frequencies depending on the task. If we drop the residual connection, W is less powerful and the only asymptotic regime is over-smoothing (LFD), independent of the spectrum of W.
5.2 Frequency response and heterophily
Negative eigenvalues flip the edge signs. Let W #ddia(() be the eigendecomposition of W yiclding the Fourier cocfficients Z(t)= F. We rewrite the discretized gradient flow F(t + r) +(t)+ (tF)) in the Fourier domain of W as Z(t (()++)= Z(t) AAt())diag(()) and note that along the eigenvectors of W, if μ,, < 0, then the dynamics is equivalent to flipping the sign of the edges. Therefore, negative edge-weight mechanisms as those proposed in Bo et al. (2021); Yan et al. (2021) to deal with heterophilic graphs, can in fact be simply achieved by a residual graph convolutional model where the mixing matrix W has negative eigenvalues. We refer to (33) in the appendix for more details.
Can other 'time-contimouss MPNNs be HFD? A framework that cannot enhance the high frequencies may struggle on heterophilic datasets. While by Theorem 5.1 we know that gradient-flow MPNNs can be HFD,
9

---

Published in Transactions on Machine Learning Research (MM/YYY)
prove that some existing time-contimuous MPNNs CGNN (Xhonneux et al., 2020), GRAND (Chamberlain et al., 2021a), PDE GCND Fliasof et al. (2021) are never HFD; in fact, in our experiments (Table 3) we validate that these models are more vulnerable to heterophily for an explicit statement including convergence rates we refer to Theorem B.4.
Theorem 5.4 (Informal). CGNN, GRAND and PDE - GCNp induce smoolhing and are never HFD.
Summary of the theoretical results. We have analysed a class of linear MPNNs that represent discrete gradient flows of Eg. This provides a 'multi-particle" interpretation for graph convolutions and sheds light onto the dynamics they generale. We have shown that the interaction between the eigenvectors and spectra of W and A is what drives the dynamics. We have also proven how such interaction ultimately leads to two opposite asymptotic regimes referred to as over-smoothing (LFD) and over-sharpening (HFD). This has allowed us to draw connections with more recent frameworks, proving that simple convolutional models are not necessarily bound to induce smoothing among features, provided that we have a residual connection, and can indeed also magnify the high frequencies, a property that may be desirable on heterophilic tasks, ultimalely providing affirmative answers to Q.1 and Q.2 in Section 1.
6 Extending the analysis to non-linear layers: the family of energy-dissipating MPNNs
As discussed above, analysing energies along MPNNs valuable approach for investigating their dynamics.
Cai & Wang (2020); Bodnar et al. (2022) showed that, under some assumptions, e is decreasing (expo- nentially) along some classes of graph convolutions, implying over-smoothing Consider the general family of comeecooninuuus) graph convolutions as in (1), with o a nonlinear activation:
also Rusch et al. (2022).
F()) (- (- F() AF( F(O)W).
(15)
Although this is no longer a gradient flow (unless a is linear), we prove that if the weights are symmetric, then the energy Eg in (6) still decreases along (15).
Thcorem 6.1. Consider o ::R. R + R satisfying x + 20(x)) 2 0. F solves (15) with Ω W being symmetric, then 1+ E (F()) is decreusing. If we discrelize the system with step size T and let e denole the most positive eigenvalue of Q & I, W 3 A if no positive eigenvalues exists take c = then
- S ((t F(112
An important consequence of Theorem 6.1 is that for non-linear graph convolutions with symmetric weights, the physics interpretation is preserved since the same multi-particle energy &g in (9) dissipates along the features. Note that in the time-discrete setting, the incquality can be interpreted as Lipschitz regularity result. In general, in the nonlinear case we are no longer able to derive exact results in the limit of many layers, yet the channel- shannel-mixing W still induces attraction attraction/repulsiion /repulsion along edges via its positive/negative cigenvalues (sce Lemma D.1). Theorem 6.1 differs from Cai & Wang (2020); Bodnar et al. (2022) in two ways: (i) It asserts monotonicity of an energy Es more general than Dir since is parametric and in fact also able to enhance the high frequencies; (ii) it holds for an infinite class of non-linear activations (beyond ReLU).
7 Experimental validation of the theoretical results
In this Section we validate the thcoretical results through ablations and real-world experiments on node classification tasks. First, we study a linear, gradient-fow MPNN of the form:
F(t 1) = F(t) + AF(t) Ws = (W + W 1/2.
(16)
where We E Rdxd is a learnable weight matrix. Thanks to Proposition 4.1 we know that this is a discretiizd gradient flow that minimizes Eg in (6). In particular, from Theorem 5.1 we derive that the positive eigenvalues of Ws induce attraction along the edges, while the negative eigenvalues generate repulsion. We note that while the uPNN-update is linear, the map associating a label to each node based on its input feature, is actually non-linear, because the encoder and decoder are, typically, non-linear MLPs. Since this family of
10

---

Published in Transactions on Machine Learning Research
Dataset
Texas
W Wiscousin
Cornell
Film
Squireel
Chameleon
Citeseer
Pubmed
Cora
Homophily
0.11
0.21
0.30
0.22
0.22
0.23
0.74
0.80
0.81
Nodes
183
251
183
7,600
5,201
2,277
3,327
18,717
2,708
#Edges
295
466
280
26,752
198, 193
31,421
1,076
41,327
5,278
GCN
60.81 + 4.4
51.57
3.51
1.75
29.96
0.84
12.77
2.38
62.63 2.02
75.9 1.4
87.31
0.41
5
80.54
4.95
74.82
33.72
0.76
51.02
1.688
68.37
1.72
88.35
87.12
0.61
GCN,/
84.86
77.3
7.85
35.13
0.61
50.84
1.92
0.42
87.79 0.93
Table 1: Comparison of GCN and the models in (16) and (17) over datasets of varying homophily.
graph convolutions amount a gradient-flow, simplified (i.e. linear) GCN with residial--onmeetion,, we adopt the notation SGCN,f (we recall that A = D 1-11/2D AD AD-1/2).
By Theorem 6.1, if we 'activate' (16) as
1) ( (AF()))s) Ws = (W + W 1/2
(17)
with a s.t. (x) 2 O, then &a in (6) is decreasing, so that we can think of such equations as more general "approximate" gradient flows, With slight abuse of notations, we refer to (17) as GCN since this is just a GCN-model with a residual connection and symmetric weights, which cnsure the monotonicity of the enorgy.
Real-world experiments: performance on node-classification tasks. By Theorem 5.1 and The- orem 6.1, both (16) and (17) can induce repulsion along the edges through the negative cigenvalues of Ws- Therefore, these models should be more robust to the heterophily of the graph when compared to the classical implementation of GCN in fact, note that linear GCN with symmetric weights is always LFD per Theorem 5.3, and that non-lincar GCN is bound to over-smooth if the singular values are sufficiently small (Oono & Suzuki, 2020; Cai & Wang, 2020). To validate this point, we run all three models on node classification tasks defined over graphs with varying homophily (Sen et al., 2008; Rozemberozki et al., 2021; Pei et al., 2020) (details in Appendix E). Training, validation and test splits are taken from Pei et al. (2020) for all datasets for comparison. The results are reported in Table 1. For all the datasets shown, we performed simple, grid search over the space m (2, 42) recall that is the depth learning rate E (0.001, and decay E {0.0005, 0.005, 0.05J. We find that both (16) and (17) surpass the performance of GCN, often by significant margin, across all heterophilic datasets. We emphasize that this is in support of our theoretical findings, since we proved that thanks to a residual connection, models enable the channel-mixing to also generate repulsion along the edges through its negative eigenvalues hence compensating for the underlying heterophily of the graphs. In particular, the results also show that imposing symmetrization over the weights has no negative impact on performance.
Remark. It is possible to significantly improve the results by considering additional Ω and W -teems in (10) and tuning, achieving strong results that often surpass more complicated benchmarks specifically designed for helerophálic lasks. Our resulls then further bring aboul questions about the qualily of these lasks and whether sraph-convolutionss models really need to be augmented or replaced for heterophálic tasks. Since the main goal of the Section amounts to validating the theory, we report these results in Appendir E.
Spectral analysis. To corrobo- rate our analysis, we also investigate spectral properties of the weight ma- trix Ws. We report statistics of the gradicnt-flow MPNN in (16) after training, over 4 datasets (two hel- erophilic ones, and two homophilic ones). We recall that fa-1 and μl0 are the most positive and most nega- tive eigenvalues of Ws, respectively.
Dataset
Squirrel
Chameleon
Pubmed
Cora
Homophily
0.22
0.23
0.8
0.81
-Ha-1/PP &D (F(O)) F(
0.68
0.71
2.43
9.58
0.83
0.65
0.35
0.19
Table 2: Spectral properties of the weights of SGCN,, and behaviour of the (normalized) Dirichlet energy.
11

---

Published in Transactions on Machine Learning Research (MM/ MYYY)
2.00
1.75
camel 60
1.50
90
1.50
10
1.25
8
2 1.00
UON ol pirt-10
0.75
0.75
0.50
0.25
0.25
0.00
0.00
5
10
20
46
50
60
tme
tim
Figure 2: Rayleigh Quotient of Ia @ A (i.e. Figure 3: Rayleigh quotient of La @ A on a com- Dir (F(t))/[[((t)l on real graphs for SGCNgf. plete bipartite graph for SGCNg5 and GCN.
In fact, by Theorem 5.1 we know that the ratio -Me-1PP indicates whether the underlying dynamics is dominated by the high frequencies (ratio smaller than 1), or by the low frequencies (ratio larger than one). We confirm this behaviour in Table 2, where we sce that the quantity -44-1/144 is larger on the homophilic graphs and is instead smaller on the heterophilic graphs, where the model generally learns to enhance the high-frequency components more. We also monitor the profile of the normalized Dirichlet energy. We see that homophilic graphs the low frequencies at the final layer is much lower than the one evaluated over the input features. This validates that our characterization of over-smoothing (LFD dynamics) and over-sharpening (HFD dynamics) is appropriate to study the dynamics of GNNs.
indeed dominant since the normalized Dirichlet energy
Over-smoothing and over-sharpening. We also validate our theoretical results that linear gradient-flows MPNNs as in (16) admit two, opposite, behaviours in the limit of many layers, depending the eigenvalues of w: (i) over-smoothing, i.e. LFD dynamics where the normalized Dirichlet energy approaches zero and (ii) over-sharpening, i.e. HFD dynamics where the normalized Dirichlet energy approaches the largest eigenvalue of the normalized Laplacian A E [1,2). Accordingly, we report the values of the normalized Dirichlet energy (i.e. the Rayleigh quotient La & A, simply denoted by RQ) over several worldworld datasets; here we are interested in validating that both over-smoothing and over-sharpening occur and that no other intermediate case is possible. Figure 2 confirms our theoretical analysis, showing that after training, the normalized Dirichlet energy approaches zero in 7 out of cases (LFD dynamics, i.e. over-smoothing) whereas it converges to its maximum, depending the graph, on the remaining two cases (HFD) setting, i.e. over-sharpening), and that no other intermediate case arises. We point out that in this experiment we study the limiting behaviour of the underlying MPNNs and not the actual performance, which will of course degrade as we increase the number of layers for both weershhothhigg and sharpenar...in. To further support our claims, we also consider the synthetic case of a complete, bipartite graph with maximal heterophily (i.c. each node is connected to all other nodes with different label). In this extreme case, pver-shappenin does actually provide the optimal classifier since the eigenvector of the graph Laplacian associated with largest frequency separates the two communities perfectly; in Figure 3 see that for (16) the dynamics is HFD (i.e. over-sharpening), while a standard GCN over-smooths due to the lack of a residual connection as per Theorem 5.3.
8 Conclusions
Summary and the messages of the work. In this work we have shown that graph-convolational models can enhance the high-frequency components of the features via the negalive eigenvalues of the weight matrices, provided that there is a residual connection. particular, we have studied a simple class of graph convolutions that represent the gradient flow of a parametric functional Eg. The dynamics of gradient-flow MPNNs is fully characterized in terms of the interactions between the eigenvalues of the (symmetric) weight matrices and those of the graph Laplacian, ultimately leading to either over-smoothing or over-sharpening, in the limit of
12

---

Published in Transactions on Machine Learning Research (MM/YYYY)
many layers. We have also extended some of the conclusions to gradient flows activated by non-lincar maps and have finally validated our analysis through ablation studies and real-world experiments.
Our work suggests that: (i) It is possible to shed light on the dynamics of GNNs once we identify an energy function that decreasing along the features, thereby providing a general framework for analyzing existing models as well as building more principled ones; (ii) graph convolutions are not simple low-pass filters but can also learn generate edge-wise repulsion once we have a residual connection; (iii) the role of heterophily as a brought into question, given that simple GCN models with residual connections already par, better than GNN models designed specifically to deal with heterophily; (iv) the idca that over-smoothing is the dominant or only asymptotic behaviour for very deep MPNNs is questioned, considering that graph convolutions can also incur over-sharpening if the negative eigenvalues of the weight matrices are sufficienlly large.
major issue for GNNs and/or the validity of benchmarks commonly adopted in node classification,
Limitations and future work. We limited our attention to a class of energy functionals whose gradient flows give rise evolution equations of the praph-convolutional type. In future work, we plan study other families of energies that generalize different GNN architectures and provide new models that more 'physics'-insse -nnspire;; we believe this to be powerful paradigm to design new GNNs. To the best of our knowledge, our analysis is first step inlo studying the interaction of the graph and 'channel-mxxing spectra; it will be interesting to explore more general dynamics (i.e., that are neither LFD nor HFD).
Acknovledgements. We thank the anonymous reviewers for valuable suggestions and feedback.
13

---

Published in Transactions on Machine Learning Research (MM/YYY)
References
Uri Alon and Eran Yahav. On the bottleneck of graph neural networks and its practical implications. In International Conference Learning Representations, 2021. URL https:/ bttta:://openrvvee...o wttpr//wwwwee/oommmiiiieee 1800Ph0CVH2.
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 Represenlations, 2020.
Lukas Biewald. Experiment tracking with weights and biases, 2020. URL http://wwwwwwwdbco... com/. Software available from wandb.com.
Marin Biloš, Johanna Sommer, Syama Sundar Rangapuram, Tim Jamuschowski, and Stephan Gün nemann. Neural flows: Efficient alternative neural odes. In Advances Neural Information Processing Systems, volume 34, 2021. URL https://proceedinga nwwripu b21191882aaaaa5f(dddddddt Paper pdf.
Deyu Bo, Xiao Wang, Chuan Shi, and Huawei Shen. Beyond low-frequency information in graph convolutional networks. In AAAI. AAAT 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. aiv preprint 0rXiv:2202.04571 2022.
Xavier Bresson and Thomas Laurent. Residual gated graph convnets. arXiv preprint arXi::1711.0755 2017.
Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXin:2104.13 194778 2021.
Joan Bruna, Wojeiech Zaremba, Arthur Szlam, and Yann LeCun. Spectral networks and locally connected networks graphs. In 2nd International Conference on Learning Representations, ICLR 2014, 2014.
Chen Cai and Yusu Wang. A note on over-smoothing for graph neural networks. arXiv:2006. 18918, 2020.
arXio preprint
Ben Chamberlain, James Rowbottom, Maria I Gorinova, Michael Bronstein, Stefan Webb, and Emanuele Rossi. Grand: Graph neural diffusion. In International Conference 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 diffussio on graphs. Advonces 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. In 9th International Conference on Learning Representations, TCLR 2021, 2021. URL https://openreview
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.
14

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Francesco Di Giovanni, Lorenzo Giusti, Federico Barbero, Giulia Luise, Pictro Lio, and Michael Bronstein. On over-squashing in message passing neural networks: The impact of width, depth, and topology. ar preprint Kix:22300 02941, 2023.
James Fells and Joseph H Sampson. Harmonic mappings of riemannian manifolds. American journal of mathematics, 16(1):109-160, 1964.
Moshe Eliasof, Eldad Haber, and Eran Treister. Pde-gen: Novel architectures for graph neural networks motivated by partial diffcrential equations. Advances in Neural Information Processing Systems, 34, 2021.
Moshe Eliasof, Lars Ruthotto, and Eran Treister. Improving graph neural networks with lcarnable propagation operators. 2023.
Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorrc Geometric. In ICLR Workshop on Represenlation Learning on Graphs and Manifolds, 2019.
Justin Gilmer, Samuel $ Schoenholz, Patrick F Riley, Oriol Vinyals, and George Dahl. Neural message passing for quantum chemistry. In International Conference on Machine Learning, pp. 1263 1272. PMLR, 2017.
Christoph Goller and Andreas Kuchler. Learning task-dependent distributed representations by backprop- agation through structure. In Proceedings of International Conference on Neural Networks (ICNN '96,, volume 1, 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 2, pp. 234-734.
LEEE, 2005.
E. 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 Vertez- Fretue-Fre 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 in deep neural networks. Computer Vision and Image Understanding, 187:102786, 2019.
Ron Kimmel, Nir Sochen, and Ravi Malladi. From high energy physics to low level vision. In International Conference on Scale-Space Theories in 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 https:// //ppr//////eereerr netf/forur/1dd-SJ4a//y
Qimai Li, Zhichao Han, and Xiao-Ming Wu. Deeper insights into graph convolutional networks for semi- supervised learning. In Proceedings of the AAAI conference on artificial intelligence, volume 32, 2018.
Derek Lim, Felix Hohne, Xiuyu Li, Sijia Linda Huang, Vaishnavi Gupta, Omkar Bhalerao, and Ser Nam Lim. Large scale learning non-homophilous graphs: New benchmarks and strong simple methods. Advances in Neural Information Processing Systems, 34:20887 20902, 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 do node classification? arXiv preprint ar Xix::2100.00 05641, 2021.
15

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Sitao Luan, Chenqing Hua, Qincheng Lu, Jiaqi Zhu, Mingde Zhao, Shuyuan Zhang, Xiao- Wen Chang, and Doina Precup. Revisiting heterophily for graph neural networks. Advances in neural in information processing systems, 35:1362 1375, 2022.
Sunil Kumar Maurya, Xin Liu, and Tsuyoshi Murata. Improving graph noural networks with simple archilecture design. arviv preprint urXiv:2105006334 2021.
Christopher Morris, Martin Ritzert, Matthias Fey, William L Hamilton, Jan Eric Lenssen, Gaurav Rattan, and Martin Grohe. Weisfeiler and leman go neural: Higher-order graph neural networks. In AAAT Conference Artificial Intelligence, pp. 4602 4609. AAAI Press, 2019.
Hoang Nt and Takanori Machara. Revisiting graph neural networks: All we have is low-pass filters. arXiv preprinl or Xiv: Xio::195.00950, 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 DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In NeurIPS. 2019.
Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang. Geom-gen: Geometric graph convolutional networks. Tn 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):
029-639, 1990.
Oleg Platonov, Denis Kuznedelev, Michael Diskin, Artem Babenko, and Liudmila Prokhorenkova. A critical look at the evaluation of gnns under helerophily: are really making progress? ar XrXie preprinl 0rXiv:2302.11614 2023.
Benedek Rozemberezki, Carl Allen, and Rik Sarkar. Multi-scale attributed node embedding. Journal of Complex Networks, (2):cnab,1,, 2021.
T Konstantin Rusch, Benjamin Chamberlain, James Rowbottom, Siddhartha Mishra, and Michael M Bronstein. Graph-coupled oscillator networks. In International Conference Machine Learning, 2022.
T Konstantin Rusch, Michael M Bronstein, and Siddhartha Mishra. A survey oversmoothing in graph neural networks. arXie preprint arXiv:2303. 211i::3333.0099,, 2023.
Franco Scarselli, Marco Gori, Ah Chung Tsoi, Markus Hagenbuchner, and Gabriele Monfardini. The graph neural network model. IEEE transactions on neural networks, 20(1):61-80, 2008
Prithviraj Sen, Galileo Namala, Mustala Bilgie, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad. Collective classification in network data. Al magazine, 29(3):93 93, 2008
Alossandro 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.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pictro Liò, and Yoshua Bengio.
Graph attention networks. In International Conference on Learning Representations, 2018. URL https:
//openrevi ew. www.:eooumidddr..aeee....
16

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Yuelin Wang, Kai Yi, Xinliang Liu, Yu Guang Wang, and Shi Jin. Aemp: Allen-cahn message passing for graph neural networks with particle phase transition. aixiv preprint or arXin:2206.05497 2022.
Felix Wu, Amauri Souza, Tianyi Zhang, Christopher Fifty, Tao Yu, and Kilian Weinberger. Simplifying graph convolutional networks. In International conference machine learning, pp- 6861-6871. PMLR, 2019.
URL http://proceedingss mlr prwps/v97 www////w///////// html.
Louis-Pascal Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. International Conference on Machine Learning, pp. 10432-10441. PMLR, 2020. URL http://pproceeding. mlr press/ v119/shonneux20a. html.
Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? In
ICLR. OpenReview.net, 2019.
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 arXiv:2102.06. 06462, 2021.
Yu-Li You, Wenyuan Xu, Allen Tannenbaum, and Mostafa Kaveh. Bchavioral analysis of anisotropic diffusion in image processing. IEEE Transactions Image Processing, 5(11):1539 1553, 1996.
Chengxi Zang and Fei Wang. Neural dynamics on complex networks. Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 892-902, 2020.
Xitong Zhang, Yixuan He, Nathan Brugnone, Michael Perlmutter, and Matthew Hirn. Magnet: A neural network for directed graphs. Advances in neural information processing systems, 1-220033-271111 2021.
Lingxiao Zhao and Leman Akoglu. Pairnorm: Tackling oversmoothing in gnns.
ar preprint
arXiu:1909 1tt.199999222.... 12223, 2019.
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 Li, Soo-Hyun Choi, and Xia Hu. Dirichlet energy constrained learning for deep graph neural networks. Advances Newral Information Processing Systems, 34:21834 21846, 218484,411846, 2021.
Jiong Zhn, 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.
17

---

Published in Transactions on Machine Learning Research (MM/YYN)
Overview of the appendix
To facilitate navigating the appendix, where we report several additional theorctical 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 our work. We also review important elementary properties of the Kronecker product of matrices that used throughout our proofs in Appendix A.2. We finally comment on the choice of the normalization (and symmetrization) of the graph Laplacian, briefly mentioning the impact of different choices.
In Appendix B.1 we derive the energy decomposition reported in (9). In Appendix B.2 we derive additional rigorous results to justify our characlerization of LFD and HFD dynamics in Definition 4.2 along with examples. We also formalize (and prove) more quantitatively Theorem 4.3 in Theorem B.3. In Appendix B.3 we report Theorem 5.4 including convergence rates and over-smoothing results. In Appendix B.4 we explore the special case of Ω = W which is equivalent to choosing A rather than as messa pressag-ppssing matrix, hence providing new arguments as lo why propagating messages using rather than A is 'more robust'. Finally in Appendix B.5 we formally derive an analogy between the continuous energy used for manifolds (images) and subsct of the parametric energies in (5).
In Appendix C we prove Theorem 5.1, Corollary 5.2, and Theorem 5.3.
In Appendix D we prove Theorem 6.1 and an extra result confirming that even in the non-lincar case the channel-mixing W still induces attraction and repulsion via its spectrum, thus magnifying the low high frequencies respectively.
In Appendix we report additional details on the evaluation along with further real-world experiments and ablation studies.
Additional notations and conventions used throughout the appendix. Any graph G is taken to be connected. We order the eigenvalues of the graph Laplacian 0= do X1 ... An-1 An-1 S 2 with associated orthonormal basis of eigenvectors { so that in particular we have AΦd = 0. Moreover, given symmetric matrix B, we generally denote the spectrum of B by spec(B) and if B is positive semi-definite (written as B 2 0), then (B) denotes the positive, smallest eigenvalue of B. Finally, if we write (()) we always take the norm to be the Frobenius one and tacitly assume that the dynamics s.t. the solution is not zero. Without losing generality, below we always take -11///2A--/22
A Motivations and important preliminaries
A.1 Discussion on continuous Dirichlet energy and harmonic maps
Starting point: a geometric parallelism. To motivate a gradient-fow approach for GNNs, we start from the continuous case. Consider a smooth map f R" + (Rd,h) with constant metric represented by H 1 0. The Dérichlet energy of f is defined by
E(fhh) = 2 [ 11V de
n Σ
be (x)
(18)
gr==i==
and measures the 'smoothness' of f. A natural approach to find minimizers of E called harmonic maps was rntrudd in Eells & Sampson (1964) and consists in studying the gradient flow of E, wherein a given map f(0) = fo is evolved according to f(t) = -V (()( These type of evolution equations have historically been the core of variational and PDE-lused image processing; particular, gradient flows of the Dirichlet energy (Kimmel al., 1997), were shown to recover the Perona-Mall Malik nonlinear diffusion Perona & Malik (1990).
18

---

Published in Transactions on Machine Learning Research (MM/YYYY)
A.2 Review of Kronecker product and properties of Laplacian kernel
Kronecker product. In this subsection we summarize few relevant. notions pertaining to the Kronecker product of matrices, that are going to be applied throughout our spectral analysis of gradient flow equations for GNNs in both the continuous and discrete time setting.
Given a matricial equation of the form
Y = AXB,
we can vectorize X and Y by stacking columns into vec( X) and veo(Y) respectively, and rewrite the previous system as
vec Y) (B
@ A) vec(X).
(19)
If and B are symmetric with spectra (A) and spec(B) respectively, then the spectrum of B @ A is given by spec(A) spec(B). Namely, if Ax = X x and By = X >"y, for x and y non-zero vectors, then 1 is an eigenvalue of B & A with eigenvector y & x:
(B BoA)yox- (1 x
(20)
One can also define the Kronecker sum of matrices A E guxn and gdxd as
&B
(21)
with spectrum spec(A + B) (AA : \ Spec(AA) A E spec((B)).
Additional details on eir and the choice of Laplacian. We recall that the classical graph Dirichlet energy Dir is defined by
Dir trace (F).).
We can use the Kronecker product to rewrite the Dirichlet energy as
e (F) = v) T (Ia Ø A)vec(F),
(22)
from which we immediately derive that & e (F) (Ia )Fve since A is symmetric and hence recover the gradient flow in (3) leading to the graph heat equation across each channel.
Before we further comment on the characterizations of LFD and HFD dynamics, we review the main choices of graph Laplacian and the associated harmonic signals (i.e. how we can characterize the kernel spaces of the given Laplacian operator). Recall that throughout the appendix we always assume that the underlying graph G is connected. The symmetrically normalized Laplacian A I- A is symmetric, positive semi-definite with harmonic space of the form (Chung & Graham, 1997)
ker( A) span 1 : 1, (1 1)T).
(23)
Therefore, if the node features F(I) satisfy F(1) GNN) with initial condition F(0), and we assume that over-smoothing occurs, meaning that (t" + 0 for - oc for each column S r d, then the only information persisting in the asymptotic regime the degree and any dependence on the input features is lost, as studied in Oono & Suzuki (2020); Cai & Wang (2020). A slightly different behaviour occurs if instead of A, consider the unnoraalizdd Laplacian = D A with kernel given by span(1,)) meaning that if Lf (t) + 0 as + do for each < < d, then any node would be embedded into single point, hence making any scparation task impossible. The same consequence applics to the random walk Laplacian ARw = I D-IA A. In particular, we note that generally a row-slochastic matrix is not symmetric then this would in fact be doubly-stochastic and the same applies to the randon-walk Laplacian (a special exception is given by the class of regular graphs). fact, general any dynamical system governed by Aww (or simply 'A) is not the gradient flow of an energy due the lack of symmetry, as further confirmed below in (24).
it was,
19

---

Published in Transactions on Machine Learning Research (MM/YYN
B Proofs and additional details of Section 4
B.1 Attraction vs repulsion: physics-inspired framework
First, we show that we can rewrite the energy in (5) using the Kronecker product formalism. It suffices to focus on the term Q-term since the other derivations are identical. We have
E( = Et. sag
On the other hand, since (Ω @ ()voc(F) IFS T, we also find
(vec(F), 2 [))ee(F)= = fi (IF = Dagfff
which shows the equivalence claimed.
Proof of Proposition 4.1. We first note that the system in (8) can be written using the Kronecker product
vec (t)) (S2 @ I.)voo(F(t))) (W & A)voo(F(t)) (W & (.))cc(F(O)..
If this the gradient flow of + (F); then we would have
v Fe
EF(() E
S2 - W@A,
(24)
which must be symmetric due to the Hessian of function being symmetric. The latter means
Ω T 2 2 I, = (WT W) A.
Therefore, for each components 3 must have
(War WBa (Wa3 Wan)) A.
If A non-diagonal, which happens as soon as there is non-trivial edge in G, the equations are satisfied only if both Ω and W are symmetric. This shows that (8) is the gradient flow of Eo and only Q and W are symmeelric.
We now rely the spectral decomposition of W to rewrite &g explicitly in terms of attractive and repulsive interactions. If we have a spectral decomposition W = (d) we can separate the positive eigenvalues from the negative ones and write
W = Vddiag(u4) T Dddiag( y = W W
the Choleski decomposition to write W+ = + and W_ TIe_ with O+ __ E dddd Equation (9) then simply follows by direct computation: namely, if we let 10-1/2 1/2) Dij, derive
Since W+ > 0, W_ 2 0,
(F) Σ (f., (2E.) - Eay (f, Wf,)
[i. W)f) Σ LL WL) Σ
Bzy ( O.f) Σ
aiy ( e O.f f
(G. W)f)) + Σ | (VF)
12
1 Σ
0 (VF)
i.j
where we have used that Σig 1 0+f1 E. 2
20

---

Published in Transactions on Machine Learning Research (MM/YYY)
B.2 Additional details on LFD and HFD characterizations
In this subsection we provide further details and justifications for Delinition 4.2. We first prove the following simple properties.
Lemma B.1. Assume we have a (continuous) process t + F(t) E IRuxd for t 2 0. The following equivalent characterizations hold:
(1) (F)(F())) - 0 for t + oo if and only if (t) - for ssrsd.
(i) E (F() + An fort oo if and only if for any sequence t2 - oo there exist subsequence f)2 and a unit limit Fo depending the subsequence such that Af > for
1 d.
Proof. (i) Given F(!) E R
vectorize it and decompose it in the orthonormal basis {e, & 0 15 SS d, 5 t 1 with {er)-= canonical basis in R and write
vec(F((t) Σ (,),,,, @ Φb G.e(t) (vec(F()),,,,
We can then use (22) to compute the Dirichlet energy as
d n-
Dir (t))
(02 = ΣΣ (t)e 2 gap (A)
ΣΣ
Σ Σ (t),
r=
where we have used that Ao = 0 and that gap( = \1 S de for all 1 > 1. Therefore
d n-1
Dir (P())) 00
=
Geelt) -0 t0,
8=1
which occurs if and only if
LL 2 ( vec ((()))
Σ
(,)(t) -
t=1
(ii) The argument here is similar. Indeed we can write Q(1) F(t)/[F(() with Q(t) a unit-norm signal. Namely, we can vectorize it and write
vec vec(Q(t)) ((t)))
Σ eneet)e @ Φ6.
(t) 1.
r.t
r.
Then (Qi((Q(t)) + An-1 if and only if
Enee(i + Am-1, 4- oc,
rg
which holds if and only if
Σ (1)-+ 11)
(t) 0, E: de < An-1-
(25)
given the unit norm constraint. This is equivalent to the Rayleigh quotient of I, & A converging to its maximal value An-1- When this occurs, for any sequence t, - oo we have that gf < 1, mcaning that we can extract a converging subsequence that due to (25) will converge to a unit eigenvector Fo of L & A satisfying (Ia & F 1 Conversely assume for contradiction that there exists sequence
21

---

Published in Transactions on Machine Learning Research (MM/YYYY)
2 - 0o such that Ed (F() < An-1 c, for some 0. Then (25) fails to be satisfied along the
sequence, meaning that no subsequence converges to unit norm eigenvector Foc. of La 09 A with associated eigenvalue Ax-1 which is contradiction to our assumption.
Before address the formulation of (reehiggll-reuenccyddonnnna dynamics, we solve explicitly the system F(t) = (t()) in Raxd, with some initial condition F(0). We can vectorize the equation and solve (é)(F(tt) = (Ia & (F))) meaning that
vec ((())
Σ
(1- Ae) G,e(0)e, de,
G,e(0) :- (vec(F(0)), @ db).
Consider any initial condition F(O) such that
Gr.ol = Σ (vve e, @ Φo)1 > 0,
which is satisfied for each vec(F(()) E R 14,, where U- is the orthogonal complement of TR & span(
Since U+ is a lower-dimensional subspace, its complement is dense. Accordingly for a.e. F(O), we find that the solution satisfies
||voc(F(0)
Σ 2 + Ole 20ap(A)0) =
(F() - () O(e (()=)))((A)))
with P) the projection onto Rd @ ker (e)(A), We see 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 (e)(F())
P.l ker (A) vec(F(()) + Ole -aa(()))(( Paeeee kerl A vve(FF(O) - vec(Fxo),
Voc (()))]
(I vec (Fc)FF O(e ((e-2xmp(A)1))
such that (I. A Feo 0 which means Ar 0, for each column 1 d. Equivalently, one can compute E EDi (F(t) and conclude that the latter quantity converges to zero t - co by the very same argument.
In fact, this motivates further the nomenclature LFD and HFD. Without loss of generality we focus now on the high-frequency case. Assume that we have a HFD dynamics i + F(!), ie. DDi (F)) - Aa-1, then we can vectorize the solution and write vec(F(t)) = (vec)(((t)))) for some time-dependent unit vector (t()) E Innd;
vec(Q((t)) ((())==En((t) @ de. (t)(()== 1. T.6
r.t
By Lemma B.1 and more explicitly (25), we derive that the coefficients (4,,,-)) associated with the cigevenctors o, & Φx.-, are dominant in the evolution hence justifying the name high-frequency dominated dynamics.
The next result provides theoretical justification for the characterization of low (high) frequency dominated dynamics in Definition 4.2.
Lemma B.2. Consider dynamical system F(t) GNNo with initial condition F(0).
ii) GNN₃ is LFD if and only if (Ia & (ee(F(t)
0 if and only if for each sequence t, - there F()))
F(Ol
crist a subsequence tjp + 00 and Foo (depending on the subsequence) s.t. [[(())]]]
- Foc satisfying
A fo 0, for each d.
22

---

there exists Published in Transactions on Machine Learning Research (MM/YYY)
(ti) GNN₃ is HFD if and only if for each sequence tj there erist subsequence tye 30 and Fo (depending on the subsequence) s.t. F(Gx)) - Foo satisfying A 1 for each 1 S S d. [F(())] )Il
Proof. (i) Since A (1) + for each 1 S S d if and only if (Ia @ A (F())F(()) + O, we conclude that the dynamics is LFD if and only if (Ia A) vec(F(() 0 due to (i) in Lemma B.1. Consider sequence tj + c0. Since (e,(F(t,))/ [(t))] is a bounded sequence we can extract a converging subsequence tye (/ + Feo) If the dynamics is LFD, then (Ia A) vee(F(t;))) + 0 and hence we [[
[F(O]
conclude that vec(F ker(Ia @ A). Conversely, assume that for any sequence t, subsequence tj and Foo such that F(t,,)
[[ + Fxo satisfying Afx = 0, for each 1 S S d. If for a contradiction we had e > 0 and 1j + xc such that EDi (F(4)/ 2 e for j large enough then by (i) in Lemma B.1 there exist < r S d, f= 0 and subsequence ty2 satisfying
K (ee((F(t()))
& cr 10) 5(c) 0())) 0, ([(ttt)))
meaning that there no subsequence of {tz } s.t. [LS)))c((R( + 0, providing a contradiction. (ii) This is equivalent to (ii) in Lemma B.1.
Remark. We note that in Lemma B.2 LFD dynamics does not necessarily mean that the normalized solution converges to the kernel of Ia @ A i.c. one in general has always to pass to subsequences. Indeed, we can consider the simple ezample t +> vec(F(()) == cos(t)e; & Φo, for some 1 S S d, which satisfies (t) = 0 for each but it is not a convergent function due to its oscillatory nature. Same argument applies to HED.
We will now show that (8) can lead to a HFD dynamics. To this end, we assume that w 0 so that (8) becomes F(t) = (F W. According to (9) the negative eigenvalues of W lead to repulsion. We show that the latter can induce HFD dynamics as per Definition 4.2. We let Pw be the orthogonal projection into the eigenspace of W & A associated with the eigenvalue p_ == |Mo (( (> - 1). We recall that Ho is the most negative cigenvalue of W, while pa-1 is the most positive cigenvalue of W. We define CHPD by:
HHFD == min(p- - 1d-1, |10 laaplappAA -m II A), gap( 440 II + W)(A--1-1))
Theorem B.3. If p_ > Ma-1, then F(1)= = W is HFD for a.e. F(0) and indeed we have
(F 25p- (An-1 Ph F(0) O (oo-2eee)))) 220,
and F()) converges to F E BBuxd such that Ar An- for 1 S r S d. If instead p_ pdd-1 then the dynamies is LFD and [F()] converges to Foo Rand such that A 0, for 1 - r d, exponentially fast.
Proof of Theorem B.3. Once compute the spectrum of W & A via (20), we write the solution recall that A = I, so we can express the cigenvalues of A in terms of the cigenvalues of A:
(ec(F(()) Σ the (1- (0) 2 de.
with Web, = u,, for OS d 1, where { is an orthonormal basis of eigenvectors in Rd. We can then calculate the Dirichlet energy along the solution as
g Dir (t))) (vec(((()) (La X A)vec(F(t))
eu (1-a)) r.f
We now consider two cases:
23

---

Published in Transactions on Machine Learning Research (MM/YYY)
If pr, > then pr Xe) S pd-1-
If Mr < 0, then pr -xe) I (oolAn-1 1) = p-; with eigenvectors 4, Φn-1 for each r s.t. Wtb, - path, without loss of generality we can assume that An-1 is simple eigenvalue for A. In particular, if p4, 0 and pr (1 - Xe) p-, then
(r - Xe) < max((uo|(An-22 1), Ia 1)}.
where Ar is the second most negative eigenvalue of W and An-2 the second largest eigenvalue of A. In particular, we write
An- = An-1 gap(A,, - A),
|0al = |ool - gap(I La + W).
(26)
From (i) and (ii) we derive that if 16r(1 - Nx) # p-: then
pr - de) p- min(p- Hd-1 M |oo (An- 1),p- 1121(An-1 1}}
min(p- 1dd-1 |nolgap(An-1I- A), gap([no|I+ W)(An-1 - 1)] -HHDD (27)
where we have used (26). Accordingly, if p_ > Hdd1, then
EDr (F())
20p-
An-1 Σ & , (0),
Σ
m:
(1-x((1-e)fte)
2pp
(An-1lI F(O)]22 + O(e ((e (ee-2teeDD))
By the same argument can factor out the dominant term and derive the following limit for + and for a.e. F(0) since Pw vec(F(()) = 0 only if vec(F(0)) belongs to a lower dimensional subspace of Rnd;
vec(F(t) vec (F()))
Pw vec(F(()) + ((E-"KED")(( PW )vec(F(O)))
Pw vec (()))
(1/WW vec (ec O ((e-2ennot)
| W ved F F())
where the latter is unit vector vec(F.o) satisfying (I, & A)voo(Foc) = Vve which completes the proof for the HFD case.
For the opposite case the proof can be adapted without efforts as explicitly derived in the proof Theorem 5.1.
We emphasize that the previous result includes Theorem 4.3 as a special case.
B.3 Comparison with continuous GNNs: details and proofs
Comparison with some continuous GNN models. In contrast with Theorem 4.3, we show that three main linearized continuous GNN models are cither smoothing or more generally L.D., The lincarized PDE-GCNp model Eliasof et al. (2021) corresponds to choosing W and Ω W =K(() some time-dependent family t + K(t) Rdxd,
-((]T(())) K(I) in (8), for
FPDE-GCN.(() = - (t) K(t).
The CGNN model Xhonneux et al. (2020) can be derived from (8) by setting Ω - 2, W WW I:
FOGNN(t) AF()) + +
Finally, in linearized GRAND Chamberlain et al. (2021a) a row-stochastic matrix A(F(0)) is learned from the encoding via an attention mechanism and we have
(GRAND(() - Rw F(t) - (())
We note that if A is not symmetric, then GRAND is not gradient flow.
24

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Theorem B.4. PDE GCND, CGNN and GRAND satisfy the following:
(1) PDE GCND D is smoothing model: E (FPDE-GON, (t)) I 0.
(1i) For a.e. F(O) it holds: CGNN never HFD and if we remove the source term, then eDir (FC(FFNN(()))) CGNN (C)NN(t)))) S -aap(A);
(tii) G is connected with self-loops, FGRAND(T) + a as + co, with μ² = mean(f" 15rs d.
By (ii) the source-free CGNN-evolution is LFD independent of S. Moreover, by (iii), over-smoothing occurs for GRAND. On the other hand, Theorem 4.3 shows that the negative eigenvalues of W can make the source-free gradient flow in (8) HFD. Experiments in Appendix E confirm that the gradient flow model outperforms CGNN and GRAND on heterophilic graphs.
We prove the following result which covers Theorem 5.4.
Proof of Theorem B.4. We structure the proof by following the mumeration in the statement.
(i) From direct computation we find
dE (F(t)) d
((vec(F(t)), (Ia A)vec(F(t))))
di = -vvvoc(((1) (K (1) (t) @ A ((ve)(((())) 0,
since KT (t) @ A 2 0. Note that we have used that (A x B) & D) = AC & BD.
(ii) We consider the dynamical system
F (C) AF(() + F(t)
We can write vec (et))) = Er,e (t) & Φe, leading to the system
(,) = (1 + + 0 osesn-l, 15,,,rrrr 1srsd
We solve explicitly the system
(t,e(t) c0,)(0) (109-22))
(
2 - if Af t 1+ X xe
)
(t,e(t) = (0)(1++)) otherwisc.
We see now that for a.e. F(O) the projection (Ia @ tx,, ( .) (vee(((()) is never the dominant term. fact, if there exists s.t. X > An-1, then X - de > 2 A for any other non-maximal graph Laplacian eigenvalue. It follows that there is no Q s.l. the normalized solution maxinizes the Rayleigh quotient of Ia & A, proving that CGNN is never HFD.
If we have no source, then the CGNN equation bocomes
F(t)= (F F())
= vc) (5i @ (-)))vec(F(t)),
using the Kronecker sum notation in (21). It follows that we can write the vectorizod solution in the basis ( 8 de)re
vcc(F(tt) c.
Σ
(,)(0) & Φo + Ole Sap(AAI
Σ
(0) & Φo
1=1==
r:
EG ( $b Φb ))
+ dft
Ole O(e-Bap(A)) BaP(A)t)
nexo
25

---

Published in Transactions on Machine Learning Research (MM/YYYY)
meaning that the dominant term is given by the lowest frequency component and in fact, if we normalize find ED (F()/[IF((I) S -a-
(iii) Finally, we consider the dynamical system induced by linear GRAND
FGRAND(T) = RR = -I ((F)
Since we have no channel-mixing, without. loss of generality we can assume that d == - one can then exlend the argument to any entry. We can use the Jordan form of A to write the solution of the GRAND dynamical system as
f(() Paliag(e'1 2. ( PF)),
for some invertible matrix P of eigenvectors, with
1 t
m
Jut -
where my are the eigenvalue multiplicities. Since by assumption 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 N such that )() > 0 for each entry (1,5)- Accordingly, we can apply Perron Theorem to derive that any eigenvalue of A has real part smaller than except the eigenvalue X with multiplicity one, associated with the Perron eigenvector 1, Accordingly, we find that each block edt decays to zero as f + do with the exception of the one eo associated with the Perron eigenvector. In particular, the projection of fo over the Perron eigenvector is just pln: with p 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 (8) where Ω = W, and comment on why we generally expect framework where the propagation governed by A to be more flexible than one with A. If = W and we suppress the source lerm i.e. W 0, the gradient flow (8) becomes
Ftt = (F)
(28)
We note that once vectorized, the solution to the dynamical system can be written as
d-1 n-
vec (ec(F(t))
Σ Σ
Aot
(0) @ dr.
P=0 t=
In particular, we immodiatcly deduce the following counterpart to Theorem 4.3:
Corollary B.5. If spec W) NR_ + 0, then (28) is HFD for a.c. F(O).
Differently from (8) the lowest frequency component is 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-frequency projection. In contrast, this can be done Ω = O, or equivalently one replaces -A with A providing a further justification in lerms of the interaction between graph spectrum and chamnel-mizing spectrum for why araph-convolutional models use the normalized adjacency rather than the Laplacian for propogating messages Kipf & Welling (2017).
B.5 Revisiting the connection with the manifold case
In (18) constant nontrivial metric in Rd leads the mixing of the feature channels. We adapt this idea. by considering symmetric positive semi-definite H = W W with w Rddd and using it generalize Dir by suitably weighting the norm of the edge gradients
26

---

Published in Transactions on Machine Learning Research (MM/YYYY)
d Σ q,r=1
W (F) ==
Σ
Σ di(F))EE
hor (VF'),
by 2
Σ (iJ)EE
[W(VF)
(29)
2
We note the analogy with (18), where the sum over 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,j) E E. We generally treat. W as learnable weights and study the gradient flow of WW:
- F(t) 2 2
1 V Dir (((E)((()))))))) AFF(()) W. FEW
(30)
We see that (30) generalizes (3). Proposition B.6. Let Pwa be the projection onto ker(W W). Equation (30) is smoothing since
(F()) ( 2tgap(W (a)mp("W)japp.a A [(O)] & Dir ((PW @ I, I,)vcc(F(())), t t
In fact F(t) + Foo s.t. 3 Φx E Rd: for each i E V we have (fso)i = V 2 + PWr
Proof of Proposition B.6. We can vectorize the gradient flow system in (30) and use the spectral characteri- zation of W W X in (20) to write the solution explicitly as
vec(F((1)= = Σ e (v-eee (0)b,,
(0 S
r.l
where (4+r = spee( ( W)) Rzo with associated basis of orthonormal eigenvectors given by (vrre. Then (Fi((F(t)) (vec(F(t)), (Ia & A Avvec(F(t)):
20(
Σ (0)e
Σ
2 (0)e 20(xrde)
=2,=0,
30,20
= g ((PWr @ [=)vec(((0))+ +
Σ
(0e 20(erree) X
EDi & L.)ve((F(())) + An- 2 2t w g)app( A F(0))
where we recall that Pher the projection onto ker ker(WT (WTW) W) and that by convention the index & is associated with the lowest graph frequency Ao 0 by assumption connected. This proves that the dynamics is in fact smoothing. By the very argument find that
vcc(F(t)) (I. & pker) peer)yec(F(())) (Pher & L,)voc(F(O))) t-> 00,
with pker the orthogl projection onto kerA the other terms decay exponentially zero. We first focus the first quantity, which we write
(Ia 2 pker vee((F())))) Σ (0) @ do.
which has matrix representation I T E Buxd with
doo Σ (0)
By (23) we deduce that the i-th row of I E RRux is the d-dimensional vector 2 on the term
de
Ax. We now focus
(pW LL))ve(F(()))
Σ (0) 0
27

---

Published in Transactions on Machine Learning Research (MM//YYY)
which has matrix representation Erip, (v,)) In particular, the i-th row is given by
Σ
c.(00( = (0).
This completes the proof of Proposition B.6.
Proposition B.6 implies that no weight matrix W in (30) 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), (0) e R into the eigenvectors of W W with positive eigenvalues shrink along the GNN and converge zero exponentially fast as integration time (depth)
increases.
Projections of the edge gradients (0) E IRd into the kernel of W W stay invariant.
If W has trivial kernel, then nodes with same degrees converge to the same representation hence incurring over-smoothing. 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 eigenvalues are positive even for equations which lead to residual GNNs when discretized (Chen et al., 2018). According to Proposition B.6, we do not expect (30) to succeed on heterophilic graphs where smoothing processes are generally harmful. To deal with heterophily, one needs negative eigenvalues to generate repulsive forces among adjacent features.
C Proofs and additional details of Section 5
We first explicitly report here the expansion of the discrete gradient flow in (10) after m layers further highlight how this not equivalent to single linear layer with message passing matrix A as for SGCN Wu et al. (2019). For simplicity suppress the source term.
F(t + +) = F(t) + 7 + W)
(F + +T) = (Ind + (-02 ② I, + W & A)) (ec(F(t))
m
(ec(F(mr))
f2 2 W A) vec(F(O))
(31)
and we see how the message passing matrix A actually enters the expansion after layers with each power 0 < k S m. This is not surprising, given that we are discretizing linear dynamical system, meaning that we are approzimating an exponential matriz.
C.1 From energy to evolution equations: exact expansion of the GNN solutions
We first address the proof of the main result.
Proof of Theorem 5.1. We consider a linear dynamical system
F(t+T) = F(f) + AFF()W,
with W symmetric. We vectorize the system and rewrite it as
vec(F(L + 7)) = (Ina + & A)vec(F(1))
which in particular leads to
vec(F(mT)) (Ind + -WW &
28

---

Published in Transactions on Machine Learning Research (MM/YYY)
We can then write explicitly the solution
vec(F( = (1 (1++T,(L -x)) c,e(0),, & oc.
We now verify that by assumption in (12) the dominant term of the solution is the projection into the cigenspace associated with the cigenvalue p- = | (An 1). The following argument follows the same structure in the proof of Theorem B.3 with the extra condition given by the step-size. First, we note that for any such that Hr > 0, we have
|1+ 1I Tp--1 +1++++00000 7,, Ae)
since required p- > Md-1 in (12). Conversely, if μ,, then
||+ TT,(111 Xx) S max [|1 -p- [1
Assume that T|Hol > 1, otherwise there is nothing to prove. Then 1 + T-- | 7| if and only if
T/mo|(2 An-1) < 2,
which is precisely the right inequality in (12). We can then argue exactly as in the proof Theorem B.3 to derive that for each index r such that ,4, and 1r * μo: then
+ T1, Xe) max{|1 + Tluil(An-1 2 1)1}
with in and An- defined in (26). We then introduce
ÓHPD == max (1x-1,11,,--11 | A), - (An-1 1)gap((mo|I+ W), lnol - 2
(32)
and conclude that
f((TT) = Σ (1+ xe)) (,(0))))),
= =1++Tp_) (co.n-100) (L0-)))0((00---0))))-)--11)1 40 (0((((++TTp)))
Σ c,((0))e(()))))) (1)))(1A)))p-
1
Tp- (con-
which completes the proof of (13).
Conversely, if < 14-1, then the projection onto the eigenspace spanned by tha-1 Φo is dominating the dynamics with exponential growth (1 T1 + 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 symmotric normalized Laplacian A, we have (0(1)= = / Ed,/2/El, which then yields (14).
We can now also address the proof of Corollary 5.2.
Proof of Corollary 5.2. Once have the node-wise expansion we can simply compute the Rayleigh quotient. of I A. We report the explicit details for the HFD case since the argument for LFD extends without relevant modifications. Using (11), we can compute the Dirichlet energy along solution of F(t+T) = F(t) (FF satislying (12) by
Dir (F(())
Σ (1 (1, Ae)) (0)(
= (1 + TP_) 22m An-1 Σ
A,,,,, (0) + +0(
1(11 (1 Tp-
1 (TTTTD 2m
Σ
(7((0)))))
r:
0 1 [++THHED
)):
2m
=(1 = pp_
(Nn-11
1+ TP-
29

---

Published in Transactions on Machine Learning Research (MM/YYYY)
where Pw is the orthogonal projector onto the eigenspace associated with the cigenvalue p- 1Hol(An-1- 1). In particular, since
vec(F(mT)) = (1+ veo vec(F(()) + O(6")),
we find that the dynamics is HFD with vec(F(t)//lvee(((()) converging to the unit projection of the initial projection by Pw provided that such projection is not zero, which satisfied for a.e. initial condition F(O).
Proof of Theorem 5.3. If we drop the residual connection and simply consider F(t + r) = (F) then
vCc(F(mT)) = (r W @ (cc(F)).)
Since G is not bipartite, the Laplacian spectral radius satisfies An-1 2. Therefore, for each pair of indices (r, E) we have the following bound:
[p,(1 - -c)1 S max Ho0h
and the inequality becomes strict é > 0, i.e. de > O. The eigenvalues pa-1 and po are attained along the eigenvectors 4 & Oe and &h & Φo respectively. Accordingly, the dominant terms of the evolution lie in the kernel of La A, meaning that for any Fo with non-zero projection in ker(Ia @ A) which is satisfied by all initial conditions except those belonging to lower dimensional subspace the dynamics is LFD. In fact, without loss of generality assume that |POl pd-1, then
veo(F(mTT) |ool pa
(-1)"Trr.,00))) @ Φ0
+ | mm (O(((m))(
Σ
(ebu @ @o)( @ @u) ) ((e(()))
Mo
with ((m) + 0 00, which completes the proof.
Gradient flow as spectral GNNs. We finally discuss (10) from the perspective of spectral GNNs in Balcilar et al. (2020). Let us assume that W O, Ω = 0. If we let A diag(()) be the eigendecomposition of the graph Laplacian and }ur} be the spectrum of W with associated orthonormal basis of cigenvectors given by {ırr), and we introduce z" (t) :V can rewrile the discretized gradient flow as
V + IR defined by (t) (f) bb)), then we
21 -VI+)) AA)) (1) = (1) (1),,, 0sr
(33)
Accordingly, for each projection into the r-th eigenvector of W, we have a spectral function in the graph frequency domain given by À 71, X). If μ4, we have a low-pass filter while if μ4, we have high-pass filter. Moreover, we see that along the eigenvectors of W, if Mr < 0 then the dynamics is equivalent to flipping the sign of the edge weights, which offers direct comparison with methods proposed in Bo et al. (2021); Yan et al. (2021) where some 'attentive" mechanism is proposed to learn negative edge weights based on feature information.
The previous equation simply follows from
2[(+ (E(++)),)), = (f.(t) + W(Af(t)):,)) v.)
= = + pr Σ (t)=
which concludes the derivation of (33).
30

---

Published in Transactions on Machine Learning Research (MM/YYYY
D Proofs and additional details of Section 6
Proof of Theorem 6.1. First check that if time is continuous, then Eg in (6) is decreasing, We use the Kronecker product formalism to rewrite the gradient VFEg(F) vector in R'd, explicitly, we get
(Ω & I, W & A)voc(F) + (W & I, vec (((()))()))) E
2
It follows then that
1 (F()))
(VFEo(F(t))) voc(F(t))
2
If we introduce the notation Z(t) (F((( then we can rewrite the derivative as
dF Z(t) ((tt)((Z(())) Σ {Z))
by assumption on σ. The discrete case follows similarly. Let us use the same notation as above so we can write +) =(t)) (Z(()) with Z(t) (FF)))(((t).
(F r)) (v w 8 A)vec (T) + 2(W & L.)vec(F(O)))
W A)vec(F(( 7)) 2 @ [.)ve((((O)))
(vo(F())) W & A) (v)c - ro (()))) 2(W & (F))vc(F(0)))
(voe(F()),(0)I)) W A) vcc(F 2 vccc (()))
((ve, W (Z
+ (((((()))(()))))))))) W Veec (F(() 2(W & L.)vee (F()))
+ (Z)))
A)o (Z(t)))
By using that Ω I, W & A is symmetric, we find that
= & (7)
w & A)vec + (W 8 (.))ec(F(O))
+
W A) (F(tt)
(o((())) (r 2(0) (F( F(t), ( I, w & -) - -FF)) s Eo(F(1)) +
+) -
where again we have used that T (Z) 2 0. This completes the proof.
To further support the principle that the effects induced by W are similar even this non-linear setting, we consider a simplified scenario.
Lemma D.1. If we choose Ω = w = diag(w) with wf 5 for O and W = 0 i.e. + F(()
solves the dynamical system
F(() = (-
with (x) 2 0, then the standard graph Dirichlet energy satisfies
Die dE(F)(F(t)
2 0.
dt
31

---

Published in Transactions on Machine Learning Research (MM/YYYY)
Proof. This again simply follows from directly computing the derivative:
Di (F(t))
1 d
(t 500
)
Σ
dd
Vd.
(s,))EE
Σ Σ Af ) (--))))))))))) SEV
= E E ieV
AF" (ool" 20.
Important consequence: The previous Lemma implies that even with non-linear activations, negative eigenvalues 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 more arbitrary choices is not immediate and we reserve this for future work.
E Additional details on experiments
E.1
Additional results on real-world tasks
General details. The graph-connolutionall models in (16) and (17) are implemented in PyTorch (Paszke et al., 2019), using PyTorch geometric Fey Lenssen (2019) and torchdiffeq (Chen et al., 2018). Hyperpa- rameters were tuned using wandb (Biewald, 2020) and random grid search. Experiments were run on AWS p2.8xlarge machines, cach with 8 Tesla V100-SXM2 GPUs.
We further improve our results the node-clessificttion tasks of by considering the following, more general, parameterizations. We consider an instance of the gradient-flor framework, termed GRAFF, of the form:
GRAFF F(t+ = Ft) + 7 (- FF (f))diag(w) + WF (F(()))
(34)
where w Rd and E I, W is a symmetric d x d-matrix shared across layers and (node-wise) encoder and decoder are MLPs. We descibe different implementations of W below. In accordance with Theorem 6.1, we also consider non-linear variant (termed GRAFFNL)
GRAFFNL: F(t+T) F(t)
F(t) + TO (-F(t)diag(w) + (F - (F(()),
(35)
with a s.t. (x) 2 O. We note that differently from (16) and (17), these models also allow for a source term controlled by B and a diagonal term acting on the state of a node, independent of the graph structure, controlled by w.
Methodology. In the experiments we rely on the following parameterizztions. We implement encoder ven and UDD as single linear layers or MLPs. We consider the structure of W that offer explicit control over its spectrun. In the (DD)-case, we consider a W E Raxd asl + ras and set W diag(w) w Due to the Gershgorin Theorem the eigenvalues of W belong to [wa - -E,l a3 model 'can' easily re-distribute mass in the spectrum of W via Ja, Ta- This generalizes the decomposition of W in Chen al. (2020) providing a justification in terms of its spectrum. For (D) we take W be diagonal.
domminnnt (DD) and diagonal (D) choices for
symmetric with zero diagonal and w e gd defined by Wa = Ga Σ, Wo
Wo Wagl.wa. Wa + E W Wooll so the
Real world experiments. In Table 3 we test GRAFF and GRAFF the same datasets with varying homophily adopted in Section 7 (Sen et al,, 2008; Rozemberczki et al., 2021; Pei et al., 2020). We use results provided in Yan et al. (2021, Table 1), which include GCNs models, GAT (Veličkovćć et al., 2018), PairNorm (Zhao& & Akoglu, 2019) and models designed for heterophily (GGCN (Yan el al., 2021), Geom-GCN (Pei el al., 2020), H2GCN (Zhu et al., 2020) and GPRGNN (Chien et al., 2021)). For Sheaf (Bodnar et al., 2022), a recent strong bascline with heterophily, we took the best performing variant (out of six) for each datasct. We include continuous baselines CGNN (Xhonneux al., 2020) and GRAND (Chamberlain al., 2021a) corroborate Theorem 5.4. Training, validation and test splits are taken from Pei et al. (2020) for all datasets
32

---

Published in Transactions on Machine Learning Research (M//YYYY)
Texas
Wiscousin
Cornell
Film
Squirrel
Chameleon
Citeseer
Pubmed
Cora
Hom level
0.11
0.21
0.30
0.22
0.22
0.23
0.74
0.80
0.81
A Nodes
183
251
183
7,800
5,201
2,277
3,327
18,717
2,708
# EEdges
295
466
280
26,752
195,493
31,421
4,676
44,327
5,278
# #Clasees
5
5
5
5
5
7
3
6
GGCN
84.86
4.55
86.86
3.29
85.68 + 6.63
37.54 + 1.56
55.17 + 1.58
71.14 + 1.84
77.14 + 1.45
89.15 + 0.37
87.95 + 1.05
GPRGNN
78.38
4.36
$2.94
4.21
80:27 8.11
34.63 1.22
31.61 + 1.24
46.58 1.71
77.13 1.67
87.54 0.38
87.95
1.18
H2GCN
84.86 7:23
87.65
4.98
82.70 4: 5.28
35.70 = 1.00
36.48 1: 1.86
60.11 2:15
77.11 1: 1.57
89.49 0.38
87.87
1.20
GCNII
7T.57
3.83
80.39
3.40
77.86
3.79
37.41 1.30
38.47 + 1.58
63.86 4: 3.04
77.33 1:48
90.15 0.43
88.37
1.25
Goom
GCN
2.72
64.51
3.68
60.54 3.67
31.50 1.15
38.15 + 0.92
60.00 2.881
0.47
85.35
1.57
PairNorm
60.27 + 4.34
48.43 6.14
58.92 + 3.1.5
27.40 + 1.24
50.44 + 2.04
62.74 2.82
73.50 1.47
87.53 + 0.44
86.79 + 1.01
GraphSAGE
6.14
5.58
75.95 5.001
34.23 0.99
41.61
58.73 1.68
1.30
0.50
86.90 1.04
GCN
55.14 5,16
51.76 = 3.06
60.511 + 5.30
27.32 = 1.10
53. 2 2.01
61.22 2.21
76,50 = 1.36
88.12 0.50
86,98 : 1.27
GAT
52.16 6.63
1.09
61.89 5.05
27.11 0.89
10.72 + 1.55
60.26 2.50
76.55 1.23
87.30 :..00110 1.10
86.33
MLP
80.81 4.75
85.29 = 3.31
81.89 6.40
36.53 = 0.70
28.77 + 1.56
16,21 + 2.99
71.02 # 1.90
75.69 2.00
87.16 +: 0.37
CGNN
71.35 4.05
74.31
7.26
66.22 1. 7.69
35.95 0.86
29.24 1. 1.09
46.80 1.66
76.91 1.81
57.70 0.49
87.10 1.35
GRAND
75.68 7:25
79.41
3.64
82.16 1. 7.09
35.62 1. 1.01
40.05 1. 1.111.5 1.50
54.67 1. 2.54
76.46
1.77
89.02
0.51
87.36
0.96
(heaf(mmxx)
85.95 11 5.51
89.41
1. 4.74
84.86 4.71
37.81 1. 1.15
56.34 1. 1.32
68.04 1.58
76.70
1.37
89.49
0.40
86.90
1.13
GRAFF
85.38
4.53
88.83
3.29
84.05
6.10
37.11
1.08
58.72 + 0.84
71.08
1.75
77.30
1.85
90.04
0.41
88.01
1.03
GRAFFNL
86.49 4.84
87.26
2.52
77.30 1. 3.24
35.96 0.95
59.01 1: 1.31
71.38 1.47
76.81
1.12
89.81
0.50
87.81 1.13
Table 3: Node-classffication results. Top three models are coloured by First, Second, Third
dataset
Texas
Wisconsin
Cornell
Film
Squirrel
Chameleon
Citeseer
Pubmed
Cora
Homophily
0.11
0.21
0.3
0.22
0.22
0.23
0.74
0.8
0.81
Depth
2
2
2
2
4
2
2
4
pd-1
2.07
1.33
1.37
3.33
2.48
2.66
4.1
9.17
2.49
Mo
-3.07
-0.82
-0.76
-1.45
-3.61
-3.72
-0.96
-3.89
-0.26
mean(((4)))
-0.02
0.02
0.02
0.06
-0.1
-0.1
0.13
0.15
0.11
eDir (O))
13171
5367
3687
7460
68811
25772
4749
5737
1226
CD
42674
6767
4932
12224
828532
398302
37168
119518
33700
Pirr (FO)))/(F(()))
1.68
1.33
1.35
1.1
1.77
1.86
0.55
0.82
0.57
g (F(m))
1.67
0.71
0.91
0.71
1.47
1.21
0.16
0.29
0.11
Test accuracy (mean)
81.08
81.96
72.43
34.16
52.48
68.99
76.43
88,56
87.4
Test accuracy (std)
4.52
4.45
7.13
0.9
1.7
1.53
0.51
0.97
Table 4: SGCN,, learned spectrum
for comparison. For the Cornell dataset we take the adjacency file from before the July 2022 update of Pei et al. (2020).
Results. GRAFF and GRAFFNL are both versions of graph convolutions with stronger 'inductive bias' given by the energy Ea 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. As noted in Section 7, the improved performance is mainly due to the presence of a residual connection at each layer, which simply enables the channel-mixing matrices to induce either attraction or repulsion, via their eigenvalues.
A word of caution though: these tasks are very sensitive to the hyperparameter tuning, making them generally not ideal candidates probe different models. While they suffice for our purposes, given the theoretical nature of our work, we believe them to be times very poor indicators of the 'heterophilic problem'.
33

---

Published in Transactions on Machine Learning Research
E.2 Details hyperparameters
Using wandb Biewald (2020) we performed a random grid search with uniform sampling of the continuous variables. We provide the hyperparameters that achieved the best results from the random grid search in Table Input dropout and dropout are the rates applicd to the encoder /dcoderr respectively with no dropout applied in the ODE block. Further hyperparameters decide the use of non non-lincarities, batch normalisation, parameter vector w and source term multiplier 8 which are specified the code.
w_style
Ir
decay
dropout
pippt_drroout hidden
time
step
size
chameleon
diag {lia
0.0050
0.0005
0.36
0.48
64
3.33
squirrel
diag dom
0.0065
0.0009
0.17
0.35
128
2.87
lexas
diag dom
0.0041
0.0354
0.33
0.39
64
0.6
0.5
wisconsin
diag
0.0029
0.0318
0.37
0.37
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
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
0.5
Table Selected hyperparameters for real-world datasets
E.3 Additional details the spectral experiments
We report additional details the spectrum of W in (16) and the normalized Dirichlet energy in Table 4.
We also report the distribution of the eigenvalues of the learned channel-mixing matrix Ws in Figure 4.
34

---

Published in Transactions on Machine Learning Research
Toxas homophily, 0.11
Wiscorsin homophily 0.21
Squárrel homophily. 0.22
2.0
6.5
0.0
1.5
0.75
6.53
0.25
0.50
0.25
1.02
2.8
05
W.rwal
Film homophiiyly 0.22
Chamelece homophily: 0.23
Comeil homophily: 0.3
102
6-
4-
42-
2-
22
0.0
1.5
1.5
-20
-15
-18
1.6
155
-0.50
-0.5
0.8
1.50
W.avab
W. W.svaaa
Clbsseer homophily, 0.74
homophlloo 0.8
Cora homophily 0.81
70
60-
12-
-03
5.
1.5
1.0
W.avab
W.erda
W evala
Figure 4: Histogram of SGCN of learned spectrum
35


Reviews:
[{"review_id": "sOrfaUzNbm", "reviewer": "Reviewer_azt4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**:\n* The paper is highly organized and very well written.\n* It seems that the general concept proposed here is novel and can also be applied to other MPNNs.\n* The paper is theoretically rigorous, as well as practically meaningful, because of the various experiments that show the benefit of the interpretation of MPNNs as gradient flows and the inspection of their energy.\n\n**Weaknesses / Questions**:\n* Right before the end of page 1, the authors say that \"Several works have gathered evidence that graph-convolutional models seem to struggle on some heterophilic\ngraphs\". Can you please expand on that point? Is it not true for all or most of heterophilic graphs? Are there examples of heterophilic graphs where basic MPNNs (e.g. GCN) do not struggle?\n* With respect to signed message passing architectures, the authors may consider adding a reference to \"Magnet: A neural network for directed graphs\"\n* With respect to repulsive and attractive forces, the authors may consider adding a reference to \"ACMP: Allen-Cahn Message Passing for Graph Neural Networks with Particle Phase Transition\"\n* With respect to learning low/high pass filters, the authors may consider adding references to \"Revisiting heterophily for graph neural networks\", \"Improving Graph Neural Networks with Learnable Propagation Operators\"\n* The authors call subsection 4.1 \"non-parametric\" and to 4.2 \"parametric\". I think that the terms used here may be confusing to some readers because usually the notion of \"parametric models\" is used for closed-form models (e.g., exp(-theta) where theta is a scalar parameter). Perhaps it would be useful to distinguish between the two?\n* In equation 5 the authors add a source term but I am not sure it is used later. Is it a beneficial term? Also, how is the source term here different than the residual connection made in \"Simple and Deep Graph Convolutional Networks\", and the source term in \"GRAND++: Graph Neural Diffusion with A Source Term\"?\n* If I understand correctly, the authors use the same weights W throughout all time steps. May I ask what is the influence on the theoretical results if W can change at each time step, and how may it influence the empirical results?\n* In figure 2, what is cornell_old ? is it just a typo?\n* Is it safe to assume there always exists some desired energy of an optimal state of an MPNN? If so, can such energy always be achieved using the proposed method?\n* In the appendix, the authors define the \"GRAFF\" model, but it is not shown later in table 3.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes an analysis and re-evaluation of basic MPNNs such as GCN via an energetic perspective. The paper shows that by setting the appropriate weights in models like GCN, they can be viewed as a gradient flow that is either smoothing (a commonly known observation), but also over-sharpening (a surprising result shown in this paper). \n\nThus, the authors suggest that despite the common belief that simple MPNNs like GCN are bound to always over-smooth, such issues can be avoided using their proposed parameterization.\n\nThe authors present several experiments that show the effectiveness of viewing MPNNs as gradient flows, resulting in improved performance compared to simple, known models, like GCN.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I have read the whole paper including the appendix. I enjoyed reading it, and I think that addressing my questions in my review will improve the paper.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ghtPeGmPxI", "reviewer": "Reviewer_tQzv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- The paper's results are relevant for the community, and challenge widespread beliefs by the GNN community.\n- The story of the paper is rather intuitive, going from gradient flow for time-continuous GNNs, to discrete-time GNNs (layers) and non-linearities.\n- The paper does a good job at providing explanations and intuitions for the technical results. \n- Empirical results corroborate the efficacy of the parametrization studied.\n\n**Weaknesses**\n-  I find confusing having Q1 and Q2 in the manuscript, as to my eyes they are intimate related and Q2 may imply Q1.\n- The paper tries to be really general, yet I feel some parts are not clear, and it would benefit a lot from simplifying the formulations and extending the main explanations of the paper. For example:\n  - The formulation in Eq. (1) is quite general, yet all results in the main paper assume a simplified version of it, unless I missed something.\n  - Point (i) at the beginning of page 7 is not immediately obvious to me.\n  - Non-trivial derivations, e.g. Eq. (6), are not written in the appendix nor main paper.\n- I feel that the assumption of the adjacency matrix A being symmetric can be easily missed/forgotten, and more emphasis should be put into it. An easy way to fix it, would be to focus on _undirected_ graphs, whose directed version is a symmetric adjacency matrices.\n\n_Disclaimer:_ I have not checked the correctness of the proofs in the appendix.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper explores the expressivity power of GNNs. Specifically, it answers whether simple message-passing GNNs can show a behaviour different from that of oversmoothing, in which the node features start to homogenize as more layers are added to the GNN, irrespectively of the initial feature values. By interpreting time-continuous GNNs are gradient flows of an energy function, the authors show that certain simple GNNs can indeed show oversharpening, an effect contrary to oversmoothing. Moreover, the authors extend their results to more complex GNNs with non-linear activation functions, and empirically validate the theoretical findings, showing that a properly parametrized simple GNN can outperform tailored oversmoothing solutions, casting doubts on commonly-adopted benchmarks and solutions in the literature.   ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I think this is good work, and that the results are of interest for the community. In its current state, I wouldn't require any major changes for recommending its acceptance.\n\nHowever, I do strongly believe that it would be extremely beneficial for the possible outreach of this work to remove extra content that is not fully utilized (like additional unused parametrizations), in pursuit of strengthening the message and focus of the manuscript. As it stands right now, I think it is easily accessible to only a small fraction of the community.   ", "broader_impact_concerns": "I don't think there are any broader impact concerns for this work.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "e26EbogUnL", "reviewer": "Reviewer_fVAu", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n\n- Based on normalized Dirichlet energy, over-smoothing and its opposite state (called over-shaping in this paper) are defined in a unified manner (Definition 4.2.)\n- In the previous analysis of over-smoothing, only the scale of the GNN (specifically, the singular value) was considered. This paper analyzes more precisely the role of the weight matrix and finds that the negative eigenvalue is helpful for avoiding over-smoothing (Theorem 4.3)\n- The effect of residual connection on over-smoothing is further clarified (Theorem 5.3).\n\n\n**Weaknesses**\n\n- This paper first introduces the general form of the model as Equations (1) or (8). Equation (1) claims its usefulness because it includes GCN, GraphSAGE, GCNII, and GIN with residual connections as special cases. However, only its special cases are considered in most analyses. Specifically, Section 4.2.2 considers the case $\\boldsymbol{\\mathrm{\\tilde{W}}}_t=0$ and $\\boldsymbol{\\mathsf{A}}=D^{-1/2}AD^{-1/2}$. Also, Section 4.3, Theorem 5.1, and Corollary 5.2 further impose the constraint $\\boldsymbol{\\Omega}_t=0$. Among the above examples above, only GCN satisfies this parametrization. The exception is Proposition 4.1, which used the general form of (8). However, since it only provides the necessary conditions for defining a GNN as a gradient flow, I think it may not answer why we should consider the ODE of the form (1) or (8).\n- The energy $\\mathcal{E}_\\theta$ is the Diriclet energy when $\\boldsymbol{\\mathrm{\\Omega}}=\\boldsymbol{\\mathrm{W}}=\\boldsymbol{\\mathrm{I}}$ (and $\\boldsymbol{\\tilde{\\mathrm{W}}}=0$) (Equation (5)). Since the model analyzed in Theorem 4.3 assumes $\\boldsymbol{\\mathrm{\\Omega}}=\\boldsymbol{\\mathrm{W}}=0$, the model is not a gradient flow of Dirichlet energy (or at least, it is not shown explicitly). However, LFD and HFD are defined using normalized Dirichlet energy. So, the model and the LFD/HFD definitions use different energies. This fact raises questions about how the theoretical analysis (e.g., Theorem 4.3) uses the fact that the proposed model is defined as a gradient flow with $\\mathcal{E}_\\theta$.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper analyzes the behavior of node features of a GNN represented by a gradient flow or its discretized version. Continuous-time GNN is defined as a gradient flow by parametric energy determined by node features (Eq. (6)). Then, this paper introduces the concepts of low-frequency-dominant (LFD) and high-frequency-dominant (HFD), which quantitatively evaluate the smoothness of the node representation using the normalized energy. This paper shows that a GNN can be both LFD and HFD depending on the spectra of weight matrices. Furthermore, this paper introduces a residual GNN as a discretized version of the continuous-time GNN and shows similar results. Also, it is shown that the discretized GNN is always LFD when there is no residual connection. The proposed model was applied to the node classification task on various heterophily datasets to verify its consistency with the theoretical analysis.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "P.5: As explained at the beginning of Section 4, the main characteristic of gradient flow is that the energy decreases monotonically along the solution. However, this paper does not use such a property of gradient flow. Instead, the normalized energy converges to the maximum possible value under certain conditions in Theorem 4.3. I wonder how the theoretical analysis utilizes the benefit of gradient flows, such as the monotonically decreasing property.\n\nP.6: The *only if* part of Proposition 4.1 does not hold mathematically strictly speaking (i.e., if Equation (8) is the gradient flow of energy Equation (6), then, $\\boldsymbol{\\Omega}$ and $\\boldsymbol{\\mathrm{W}}$ are symmetric.) For example, in the extreme case where the graph has no edges (i.e., $\\boldsymbol{\\mathsf{A}}=0$), no constraint is imposed on $\\boldsymbol{\\mathrm{W}}$, which does not have to be symmetric in particular. Even if the only-if part of this proposition is not satisfied, it does not cause critical problems in constructing the proposed model. However, a more elaborate discussion is needed to make it a mathematically correct statement.\n\nP.8:\n> [...] (ii) graph convolutions can avoid over-smoothing through the negative eigenvalues of the channel-mixing by incurring the opposite behaviour of over-shaping [...]\n\nThis sentence implicitly assumes that $\\boldsymbol{\\mathrm{W}}$ has at least one negative eigenvalue when (12) holds. It is certainly true (if $\\mu_0\\geq 0$, then $\\lambda_{n-1}\\leq 2$ implies $\\mu_{d-1}/(\\lambda_{n-1}-1)\\geq \\mu_{d-1} \\geq \\mu_0$, which contradicts to (12).) However since it is not obvious, I suggest writing the proof explicitly.\n\nP.10, Theorem 6.1: [...] with $\\boldsymbol{\\Omega}, \\boldsymbol{W}$ symmetric, [...] -> [...] with $\\boldsymbol{\\Omega}$ and $\\boldsymbol{W}$ being symmetric, [...]\n\nP.10, Theorem 6.1: [...] let $c$ denote the most positive eigenvalue of $\\boldsymbol{\\Omega} \\otimes \\boldsymbol{I}_n - \\boldsymbol{W} \\otimes \\boldsymbol{\\mathsf{A}}$ [...]: Can we assume that there always exists positive eigenvalues?\n\nP.11:\n> We emphasize that this is in support of our theoretical findings, since we proved that thanks to a residual connection, graph-convolutional models enable the channel-mixing to also generate repulsion along the edges through its negative eigenvalues hence compensating for the underlying heterophily of the graphs.\n\nThis sentence implicitly assumes that the learned weights have negative eigenvalues. However, I wonder whether it is empirically confirmed in numerical experiments.\n\nP.11: Table 2: I think it is easier to read if we reverse the denominator and numerator of $\\frac{\\mathcal{E}^{\\mathrm{Dir}}(\\boldsymbol{\\mathrm{F}}(0))/\\|\\boldsymbol{\\mathrm{F}}(0)\\|^2}{\\mathcal{E}^{\\mathrm{Dir}}(\\boldsymbol{\\mathrm{F}}(m))/\\|\\boldsymbol{\\mathrm{F}}(m)\\|^2}$\n\nP.12: enhanche -> enhance", "broader_impact_concerns": "There are no major concerns about broader impacts.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: All 3 reviewers agree that the claims are supported by clear evidence, that the work is relevant to the TMLR community, and recommended to accept the manuscript. All the concerns and clarifications asked by the reviewers were addressed during the rebuttal and revisions. I thus recommend to accept the paper as is. 
Decision: Accept (poster)