#  $p$-Laplacian Based Graph Neural Networks

Guoji Fu  Peilin Zhao  Yatao Bian

###### Abstract

Graph neural networks (GNNs) have demonstrated superior performance for semi-
supervised node classification on graphs, as a result of their ability to
exploit node features and topological information. However, most GNNs
implicitly assume that the labels of nodes and their neighbors in a graph are
the same or consistent, which does not hold in heterophilic graphs, where the
labels of linked nodes are likely to differ. Moreover, when the topology is
non-informative for label prediction, ordinary GNNs may work significantly
worse than simply applying multi-layer perceptrons (MLPs) on each node. To
tackle the above problem, we propose a new $p$-Laplacian based GNN model,
termed as pGNN, whose message passing mechanism is derived from a discrete
regularization framework and can be theoretically explained as an
approximation of a polynomial graph filter defined on the spectral domain of
$p$-Laplacians. The spectral analysis shows that the new message passing
mechanism works as low-high-pass filters, thus rendering pGNNs effective on
both homophilic and heterophilic graphs. Empirical studies on real-world and
synthetic datasets validate our findings and demonstrate that pGNNs
significantly outperform several state-of-the-art GNN architectures on
heterophilic benchmarks while achieving competitive performance on homophilic
benchmarks. Moreover, pGNNs can adaptively learn aggregation weights and are
robust to noisy edges.

Graph Neural Networks, $p$-Laplacian

.tocmtchapter mtchaptersubsection mtappendixnone

##  1 Introduction

In this paper, we explore the usage of graph neural networks (GNNs) for semi-
supervised node classification on graphs, especially when the graphs admit
strong heterophily or noisy edges. Semi-supervised learning on graphs is
ubiquitous in a lot of real-world scenarios, such as user classification in
social media (Kipf & Welling, 2017), protein classification in biology
(Velickovic et al., 2018), molecular property prediction in chemistry
(Duvenaud et al., 2015), and many others (Marcheggiani & Titov, 2017; Satorras
& Estrach, 2018). Recently, GNNs have become the de facto choice for
processing graph structured data. They can exploit the node features and the
graph topology by propagating and transforming the features over the topology
in each layer and thereby learn refined node representations. A series of GNN
architectures have been proposed, including graph convolutional networks
(Bruna et al., 2014; Henaff et al., 2015; Niepert et al., 2016; Kipf &
Welling, 2017; Wu et al., 2019), graph attention networks (Velickovic et al.,
2018; Thekumparampil et al., 2018), and other representatives (Hamilton et
al., 2017; Xu et al., 2018; Pei et al., 2020).

Most of the existing GNN architectures work under the homophily assumption,
i.e. the labels of nodes and their neighbors in a graph are the same or
consistent, which is also commonly used in graph clustering (Bach & Jordan,
2004; von Luxburg, 2007; Liu & Han, 2013) and semi-supervised learning on
graphs (Belkin et al., 2004; Hein, 2006; Nadler et al., 2009). However, recent
studies (Zhu et al., 2020, 2021; Chien et al., 2021) show that in contrast to
their success on homophilic graphs, most GNNs fail to work well on
heterophilic graphs, in which linked nodes are more likely to have distinct
labels. Moreover, GNNs could even fail on graphs where their topology is not
helpful for label prediction. In these cases, propagating and transforming
node features over the graph topology could lead to worse performance than
simply applying multi-layer perceptrons (MLPs) on each of the nodes
independently. Several recent works were proposed to deal with the heterophily
issues of GNNs. Zhu et al. (2020) finds that heuristically combining ego-,
neighbor, and higher-order embeddings improves GNN performance on heterophilic
graphs. Zhu et al. (2021) uses a compatibility matrix to model the graph
homophily or heterphily level. Chien et al. (2021) incorporates the
generalized PageRank algorithm with graph convolutions so as to jointly
optimize node feature and topological information extraction for both
homophilic and heterophilic graphs. However, the problem of GNNs on graphs
with non-informative topologies (or noisy edges) remains open.

Unlike previous works, we tackle the above issues of GNNs by proposing the
discrete $p$-Laplacian based message passing scheme, termed as $p$-Laplacian
message passing. It is derived from a discrete regularization framework and is
theoretically verified as an approximation of a polynomial graph filter
defined on the spectral domain of the $p$-Laplacian. Spectral analysis of
$p$-Laplacian message passing shows that it works as low-high-pass
filters111Note that if the low frequencies and high frequencies dominate the
middle frequencies (the frequencies around the cutoff frequency), we say that
the filter works as a low-high-pass filter. and thus is applicable to both
homophilic and heterophilic graphs. Moreover, when $p\neq 2$, our theoretical
results indicate that it can adaptively learn aggregation weights in terms of
the variation of node embeddings on edges (measured by the graph gradient
(Amghibech, 2003; Zhou & Schölkopf, 2005; Luo et al., 2010)), and work as low-
pass or low-high-pass filters on a node according to the local variation of
node embeddings around the node (measured by the norm of graph gradients).

Based on $p$-Laplacian message passing, we propose a new GNN architecture,
called pGNN, to enable GNNs to work with heterophilic graphs and graphs with
non-informative topologies. Several existing GNN architectures, including SGC
(Wu et al., 2019), APPNP (Klicpera et al., 2019) and GPRGNN (Chien et al.,
2021), can be shown to be analogical to pGNN with $p=2$. Our empirical studies
on real-world benchmark datasets (homophilic and heterophilic datasets) and
synthetic datasets (cSBM (Deshpande et al., 2018)) demonstrate that pGNNs
obtain the best performance on heterophilic graphs and competitive performance
on homophilic graphs against state-of-the-art GNNs. Moreover, experimental
results on graphs with different levels of noisy edges show that pGNNs work
much more robustly than GNN baselines and even as well as MLPs on graphs with
completely random edges. Additional experiments (reported in Sec. F.5)
illustrate that intergrating pGNNs with existing GNN architectures (i.e. GCN
(Kipf & Welling, 2017), JKNet (Xu et al., 2018)) can significantly improve
their performance on heterophilic graphs. In conclusion, our contributions can
be summarized as below:

(1) New methodologies. We propose $p$-Laplacian message passing and pGNN to
adapt GNNs to heterophilic graphs and graphs where the topology is non-
informative for label prediction. (2) Superior performance. We empirically
demonstrate that pGNNs is superior on heterophilic graphs and competitive on
homophilic graphs against state-of-the-art GNNs. Moreover, pGNNs work robustly
on graphs with noisy edges or non-informative topologies. (3) Theoretical
justification. We theoretically demonstrate that $p$-Laplacian message passing
works as low-high-pass filters and the message passing iteration is guaranteed
to converge with proper settings. (4) New paradigm of designing GNN
architectures. We bridge the gap between discrete regularization framework and
GNNs, which could further inspire researchers to develop new graph
convolutions or message passing schemes using other regularization techniques
with explicit assumptions on graphs. Due to space limit, we defer the
discussions on related work and future work and all proofs to the Appendix.
Code available at https://github.com/guoji-fu/pGNNs.

##  2 Preliminaries and Background

Notation. Let ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$ be an
undirected graph, where ${\mathcal{V}}=\\{1,2,\dots,N\\}$ is the set of nodes,
${\mathcal{E}}\subseteq{\mathcal{V}}\times{\mathcal{V}}$ is the set of edges,
${\mathbf{W}}\in\mathbb{R}^{N\times N}$ is the adjacency matrix and
${W}_{i,j}={W}_{j,i}$, ${W}_{i,j}>0$ for $[i,j]\in{\mathcal{E}}$,
${W}_{i,j}=0$, otherwise. ${\mathcal{N}}_{i}=\\{j\\}_{[i,j]\in{\mathcal{E}}}$
denotes the set of neighbors of node $i$, ${\mathbf{D}}\in\mathbb{R}^{N\times
N}=\mathrm{diag}({D}_{1,1},\dots D_{N,N})$ denotes the diagonal degree matrix
with ${D}_{i,i}=\sum_{j=1}^{N}{W}_{i,j}$, for $i=1,\dots,N$.
$f:{\mathcal{V}}\rightarrow\mathbb{R}$ and
$g:{\mathcal{E}}\rightarrow\mathbb{R}$ are functions defined on the vertices
and edges of ${\mathcal{G}}$, respectively. ${\mathcal{F}}_{\mathcal{V}}$
denotes the Hilbert space of functions endowed with the inner product $\langle
f,\tilde{f}\rangle_{{\mathcal{F}}_{\mathcal{V}}}:=\sum_{i\in{\mathcal{V}}}f(i)\tilde{f}(i)$.
Similarly define ${\mathcal{F}}_{\mathcal{E}}$. We also denote by
$[K]=\\{1,2,\dots,K\\},\forall K\in{\mathbb{N}}$ and we use
$\|{\mathbf{x}}\|=\|{\mathbf{x}}\|_{2}=(\sum_{i=1}^{d}{x}_{i}^{2})^{1/2},\forall{\mathbf{x}}\in\mathbb{R}^{d}$
to denote the Frobenius norm of a vector.

Problem Formulation. Given a graph
${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$, each node
$i\in{\mathcal{V}}$ has a feature vector ${\mathbf{X}}_{i,:}$ which is the
$i$-th row of ${\mathbf{X}}$ and a subset of nodes in ${\mathcal{G}}$ have
labels from a label set ${\mathcal{L}}=\\{1,\dots,L\\}$. The goal of GNNs for
semi-supervised node classification on ${\mathcal{G}}$ is to train a GNN
${\mathcal{M}}:{\mathcal{V}}\rightarrow{\mathcal{L}}$ and predict the labels
of unlabeled nodes.

Homophily and Heterophily. The homophily or heterophily of a graph is used to
describe the relation of labels between linked nodes in the graphs. The level
of homophily of a graph can be measured by
${\mathcal{H}}({\mathcal{G}})=\mathbb{E}_{i\in{\mathcal{V}}}\left[\left|\\{j\\}_{j\in{\mathcal{N}}_{i},y_{i}=y_{j}}\right|\left/\left|{\mathcal{N}}_{i}\right|\right.\right]$
(Pei et al., 2020; Chien et al., 2021), where
$\left|\\{j\\}_{j\in{\mathcal{N}}_{i},y_{i}=y_{j}}\right|$ denotes the number
of neighbors of $i\in V$ that share the same label as $i$ (i.e. $y_{i}=y_{j}$)
and ${\mathcal{H}}({\mathcal{G}})\rightarrow 1$ corresponds to strong
homophily while ${\mathcal{H}}({\mathcal{G}})\rightarrow 0$ indicates strong
heterophily. We say that a graph is a homophilic (heterophilic) graph if it
has strong homophily (heterophily).

Graph Gradient. The graph gradient of an edge $[i,j],i,j\in{\mathcal{V}}$ is
defined to be a measurement of the variation of a function
$f:{\mathcal{V}}\rightarrow\mathbb{R}$222$f$ can be a vector function:
$f:{\mathcal{V}}\rightarrow\mathbb{R}^{c}$ for some $c\in{\mathbb{N}}$ and
here we use $f:{\mathcal{V}}\rightarrow\mathbb{R}$ for better illustration. on
the edge $[i,j]$.

######  Definition 1 (Graph Gradient).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ and a function
$f:{\mathcal{V}}\rightarrow\mathbb{R}$, the graph gradient is an operator
$\nabla:{\mathcal{F}}_{\mathcal{V}}\rightarrow{\mathcal{F}}_{\mathcal{E}}$
defined as for all $[i,j]\in{\mathcal{E}}$,

$(\nabla
f)([i,j]):=\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}f(j)-\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}f(i).$
(1)

For $[i,j]\notin{\mathcal{E}}$, $(\nabla f)([i,j]):=0$. The graph gradient of
a function $f$ at a vertex $i,i\in[N]$ is defined to be $\nabla f(i):=((\nabla
f)([i,1]),\dots,(\nabla f)([i,N]))$ and its Frobenius norm is given by
$\|\nabla f(i)\|_{2}:=(\sum_{j=1}^{N}(\nabla f)^{2}([i,j]))^{1/2}$, which
measures the variation of $f$ around node $i$. We measure the variation of $f$
over the whole graph ${\mathcal{G}}$ by $\mathcal{S}_{p}(f)$ where it is
defined as for $p\geq 1$,

$\displaystyle\mathcal{S}_{p}(f):={}$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\left\|(\nabla
f)([i,j])\right\|^{p}$ $\displaystyle={}$
$\displaystyle\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}f(j)-\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}f(i)\right\|^{p}.$
(2)

Note that the definition of ${\mathcal{S}}_{p}$ here is different with the
$p$-Dirichlet form in Zhou & Schölkopf (2005).

Graph Divergence. The graph divergence is defined to be the adjoint of the
graph gradient:

######  Definition 2 (Graph Divergence).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$, and functions
$f:{\mathcal{V}}\rightarrow\mathbb{R}$,
$g:{\mathcal{E}}\rightarrow\mathbb{R}$, the graph divergence is an operator
$\mathrm{div}:{\mathcal{F}}_{\mathcal{E}}\rightarrow{\mathcal{F}}_{\mathcal{V}}$
which satisfies

$\langle\nabla f,g\rangle=\langle f,-\mathrm{div}g\rangle.$ (3)

The graph divergence can be computed by

$(\mathrm{div}g)(i)=\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left(g([i,j])-g([j,i])\right).$
(4)

Fig. 4 in Sec. E.1 gives a tiny example of illustration of graph gradient and
graph divergence.

Graph $p$-Laplacian Operator. By the definitions of graph gradient and graph
divergence, we reach the definition of graph $p$-Laplacian operator as below.

######  Definition 3 (Graph $p$-Laplacian333Note that the definition adopted
is slightly different with the one used in Zhou & Schölkopf (2005) where
$\|\cdot\|^{p-2}$ is not element-wise and the one used in some literature such
as Amghibech (2003); Bühler & Hein (2009), where
$(\Delta_{p}f)(i)=\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left|f(i)-f(j)\right|^{p-2}\left(f(i)-f(j)\right)$
for $p>1$ and $p=1$ is not allowed.).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ and a function
$f:{\mathcal{V}}\rightarrow\mathbb{R}$, the graph $p$-Laplacian is an operator
$\Delta_{p}:{\mathcal{F}}_{\mathcal{V}}\rightarrow{\mathcal{F}}_{\mathcal{V}}$
defined by

$\Delta_{p}f:=-\frac{1}{2}\mathrm{div}(\|\nabla f\|^{p-2}\nabla f),\text{ for
}p\geq 1.$ (5)

where $\|\cdot\|^{p-2}$ is element-wise, i.e. $\|\nabla
f(i)\|^{p-2}=(\|(\nabla f)([i,1])\|^{p-2},\dots,\|(\nabla f)([i,N])\|^{p-2})$.

Substituting Eq. 1 and Eq. 4 into Eq. 5, we obtain

$\displaystyle(\Delta_{p}f)(i)={}$
$\displaystyle\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\|(\nabla
f)([j,i])\|^{p-2}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}f(i)\right.$
$\displaystyle\qquad\left.-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}f(j)\right).$ (6)

The graph $p$-Laplacian is semi-definite: $\langle
f,\Delta_{p}f\rangle=\mathcal{S}_{p}(f)\geq 0$ and we have

$\left.\frac{\partial\mathcal{S}_{p}(f)}{\partial
f}\right|_{i}=p(\Delta_{p}f)(i).$ (7)

When $p=2$, $\Delta_{2}$ is refered as the ordinary Laplacian operator and
$\Delta_{2}={\mathbf{I}}-{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}$
and when $p=1$, $\Delta_{1}$ is refered as the Curvature operator and
$\Delta_{1}f:=-\frac{1}{2}\mathrm{div}(\|\nabla f\|^{-1}\nabla f)$. Note that
Laplacian $\Delta_{2}$ is a linear operator, while in general for $p\neq 2$,
$p$-Laplacian is nonlinear since $\Delta_{p}(af)\neq a\Delta_{p}(f)$ for
$a\in\mathbb{R}$.

##  3 $p$-Laplacian based Graph Neural Networks

In this section, we derive the $p$-Laplacian message passing scheme from a
$p$-Laplacian regularization framework and present pGNN, a new GNN
architecture developed upon the new message passing scheme. We demonstrate
that $p$-Laplacian message passing scheme is guaranteed to converge with
proper settings and provide an upper-bounding risk of pGNNs in Sec. C.1.

###  3.1 $p$-Laplacian Regularization Framework

Given an undirected graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ and a
signal function with $c$ ($c\in{\mathbb{N}}$) channels
$f:{\mathcal{V}}\rightarrow\mathbb{R}^{c}$, let
${\mathbf{X}}=({\mathbf{X}}_{1,:}^{\top},\dots,{\mathbf{X}}_{N,:}^{\top})^{\top}\in\mathbb{R}^{N\times
c}$ with ${\mathbf{X}}_{i,:}\in\mathbb{R}^{1\times c},i\in[N]$ denoting the
node features of ${\mathcal{G}}$ and
${\mathbf{F}}=({\mathbf{F}}_{1,:}^{\top},\dots,{\mathbf{F}}_{N,:}^{\top})^{\top}\in\mathbb{R}^{N\times
c}$ be a matrix whose $i^{th}$ row vector
${\mathbf{F}}_{i,:}\in\mathbb{R}^{1\times c},i\in[N]$ represents the function
value of $f$ at the $i$-th vertex in ${\mathcal{G}}$. We present a
$p$-Laplacian regularization problem whose cost function is defined to be

${\mathcal{L}}_{p}({\mathbf{F}}):=\min_{{\mathbf{F}}}\;\mathcal{S}_{p}({\mathbf{F}})+\mu\sum_{i=1}^{N}\|{\mathbf{F}}_{i,:}-{\mathbf{X}}_{i,:}\|^{2},$
(8)

where $\mu\in(0,\infty)$. The first term of the right-hand side in Eq. 8 is a
measurement of variation of the signal over the graph based on $p$-Laplacian.
As we will discuss later, different choices of $p$ result in different
smoothness constraint on the signals. The second term is the constraint that
the optimal signals
${\mathbf{F}}^{*}=\operatorname*{arg\,min}{\mathcal{L}}_{p}({\mathbf{F}})$
should not change too much from the input signal ${\mathbf{X}}$, and $\mu$
provides a trade-off between these two constraints.

Regularization with $p=2$. If $p=2$, the solution of Eq. 8 satisfies
$\Delta_{2}{\mathbf{F}}^{*}+\mu({\mathbf{F}}^{*}-{\mathbf{X}})=\mathbf{0}$ and
we can obtain the closed form (Zhou & Schölkopf, 2005)

${\mathbf{F}}^{*}=\mu(\Delta_{2}+\mu{\mathbf{I}}_{N})^{-1}{\mathbf{X}}.$ (9)

Then, we could use the following iteration algorithm to get an approximation
of Eq. 9:

${\mathbf{F}}^{(k+1)}=\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}{\mathbf{F}}^{(k)}+\beta{\mathbf{X}},$
(10)

where $k$ represents the iteration index, $\alpha=\frac{1}{1+\mu}$ and
$\beta=\frac{\mu}{1+\mu}=1-\alpha$. The iteration converges to a closed-form
solution as $k$ goes to infinity (Zhou et al., 2003; Zhou & Schölkopf, 2005).
We could relate the the result here with the personalized PageRank (PPR) (Page
et al., 1999; Klicpera et al., 2019) algorithm:

######  Theorem 1 (Relation to personalized PageRank (Klicpera et al., 2019)).

$\mu(\Delta_{2}+\mu{\mathbf{I}}_{N})^{-1}$ in the closed-form solution of Eq.
9 is equivalent to the personalized PageRank matrix.

Regularization with $p>1$. For $p>1$, the solution of Eq. 8 satisfies
$p\Delta_{p}{\mathbf{F}}^{*}+2\mu({\mathbf{F}}^{*}-{\mathbf{X}})=\mathbf{0}$.
By Eq. 6 we have that, for all $i\in[N]$,

$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}}}\|(\nabla
f^{*})([j,i])\|^{p-2}\left(\frac{1}{\sqrt{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{*}-\frac{1}{\sqrt{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{*}\right)$
$\displaystyle\qquad\qquad\qquad\qquad+\frac{2\mu}{p}\left({\mathbf{F}}_{i,:}^{*}-{\mathbf{X}}_{i,:}\right)=\mathbf{0}.$

Let ${\mathbf{M}}^{(k)}\in\mathbb{R}^{N\times N}$,
$\bm{\alpha}^{(k)}=\mathrm{diag}(\alpha_{1,1}^{(k)},\dots,\alpha_{N,N}^{(k)})$,
$\bm{\beta}^{(k)}=\mathrm{diag}(\beta_{1,1}^{(k)},\dots,\beta_{N,N}^{(k)})$.
Based on which we can construct a similar iterative algorithm to obtain a
solution (Zhou & Schölkopf, 2005): for all $i\in[N]$,

${\mathbf{F}}_{i,:}^{(k+1)}=\alpha_{i,i}^{(k)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\beta_{i,i}^{(k)}{\mathbf{X}}_{i,:},$
(11)

for all $i,j\in[N]$,

${M}_{i,j}^{(k)}={W}_{i,j}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|^{p-2},$
(12)
$\alpha_{i,i}^{(k)}=1\big{/}\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+\frac{2\mu}{p}\right),\quad\beta_{i,i}^{(k)}=\frac{2\mu}{p}\alpha_{i,i}^{(k)}.$
(13)

Note that when
$\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|=0$,
we set $M_{i,j}^{(k)}=0$. It is easy to see that Eq. 10 is the special cases
of Eq. 14 with $p=2$.

######  Remark 1 (Discussion on $p=1$).

For $p=1$, when $f$ is a real-valued function ($c=1$), $\Delta_{1}f$ is a step
function, which could make the stationary condition of the objective function
Eq. 8 become problematic. Additionally, $\Delta_{1}f$ is not continuous at
$\|(\nabla f)([i,j])\|=0$. Therefore, $p=1$ is not allowed when $f$ is a real
value function. On the other hand, note that there is a Frobenius norm in
$\Delta_{p}f$. When $f$ is a vector-valued function ($c>1$), the step function
in $\Delta_{1}f$ only exists on the axes. The stationary condition will be
fine if the node embeddings ${\mathbf{F}}$ are not a matrix of vectors that
has only one non-zero element, which is true for many graphs. $p=1$ may work
for these graphs. Overall, we suggest to use $p>1$ in practice but $p=1$ may
work for graphs with multiple channel signals as well. We conduct experiments
for $p>1$ ($p=1.5,2,2.5$) and $p=1$ in Sec. 5.

###  3.2 $p$-Laplacian Message Passing and pGNN Architecture

$p$-Laplacian Message Passing. Rewrite Eq. 11 in a matrix form we obtain

${\mathbf{F}}^{(k+1)}=\bm{\alpha}^{(k)}{\mathbf{D}}^{-1/2}{\mathbf{M}}^{(k)}{\mathbf{D}}^{-1/2}{\mathbf{F}}^{(k)}+\bm{\beta}^{(k)}{\mathbf{X}}.$
(14)

Eq. 14 provides a new message passing mechanism, named $p$-Laplacian message
passing.

######  Remark 2.

$\bm{\alpha}{\mathbf{D}}^{-1/2}{\mathbf{M}}{\mathbf{D}}^{-1/2}$ in Eq. 14 can
be regarded as the learned aggregation weights at each step for message
passing. It suggests that $p$-Laplacian message passing could adaptively tune
the aggregation weights during the course of learning, which will be
demonstrated theoretically and empirically in the later. $\beta{\mathbf{X}}$
in Eq. 14 can be regarded as a residual unit, which helps the model escape
from the oversmoothing issue (Chien et al., 2021).

We present the following theorem to show the shrinking property of
$p$-Laplacian message passing.

######  Theorem 2 (Shrinking Property of $p$-Laplacian Message Passing).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$ with
node features ${\mathbf{X}}$,
$\bm{\beta}^{(k)},{\mathbf{F}}^{(k)},{\mathbf{M}}^{(k)},\bm{\alpha}^{(k)}$ are
updated accordingly to Equations 11, 12 and 13 for $k=0,1,\dots,K$ and
${\mathbf{F}}^{(0)}={\mathbf{X}}$. Then there exist some positive real value
$\mu>0$ which depends on ${\mathbf{X}},{\mathcal{G}},p$ and $p>1$ such that

${\mathcal{L}}_{p}({\mathbf{F}}^{(k+1)})\leq{\mathcal{L}}_{p}({\mathbf{F}}^{(k)}).$

Proof see Sec. D.2. Theorem 2 shows that with some proper positive real value
$\mu$ and $p>1$, the loss of the objective function Eq. 8 is guaranteed to
decline after taking one step $p$-Laplacian message passing. Theorem 2 also
demonstrates that the iteration Equations 11, 12 and 13 is guaranteed to
converge for $p>1$ with some proper $\mu$ which is chosen depends on the input
graph and $p$.

pGNN Architecture. We design the architecture of pGNNs using $p$-Laplacian
message passing. Given node features ${\mathbf{X}}\in\mathbb{R}^{N\times c}$,
the number of node labels $L$, the number of hidden units $h$, the maximum
number of iterations $K$, and ${\mathbf{M}}$, $\bm{\alpha}$, and $\bm{\beta}$
updated by Equations 12 and 13 respectively, we give the pGNN architecture as
following:

$\displaystyle{\mathbf{F}}^{(0)}=\text{ReLU}({\mathbf{X}}\Theta^{(1)}),$ (15)
$\displaystyle{\mathbf{F}}^{(k+1)}=\bm{\alpha}^{(k)}{\mathbf{D}}^{-1/2}{\mathbf{M}}^{(k)}{\mathbf{D}}^{-1/2}{\mathbf{F}}^{(k)}+\bm{\beta}^{(k)}{\mathbf{F}}^{(0)},$
(16)
$\displaystyle{\mathbf{Z}}=\text{softmax}({\mathbf{F}}^{(K)}\Theta^{(2)}),$
(17)

where $k=0,1,\dots,K-1$, ${\mathbf{Z}}\in\mathbb{R}^{N\times L}$ is the output
propbability matrix with ${Z}_{i,j}$ is the estimated probability that the
label at node $i\in[N]$ is $j\in[L]$ given features ${\mathbf{X}}$ and graph
${\mathcal{G}}$, $\Theta^{(1)}\in\mathbb{R}^{c\times h}$ and
$\Theta^{(2)}\in\mathbb{R}^{h\times L}$ are the first- and the second-layer
parameters of the neural network.

######  Remark 3 (Connection to existing GNN variants).

The message passing scheme of pGNNs is different from that of several GNN
variants (say, GCN, GAT, and GraphSage), which repeatedly stack message
passing layers. In contrast, it is similar with SGC (Wu et al., 2019), APPNP
(Klicpera et al., 2019), and GPRGNN (Chien et al., 2021). SGC is an
approximation to the closed-form in Eq. 9 (Fu et al., 2020). By Theorem 1, it
is easy to see that APPNP, which uses PPR to propagate the node embeddings, is
analogical to pGNN with $p=2$, termed as 2.0GNN. APPNP and 2.0GNN work
analogically and effectively on homophilic graphs. 2.0GNN can also work
effectively on heterophilic graphs by letting $\Theta^{(2)}$ be negative.
However, APPNP fails on heterophilic graphs as its PPR weights are fixed
(Chien et al., 2021). Unlike APPNP, GPRGNN, which adaptively learn the
generalized PageRank (GPR) weights, works similarly to 2.0GNN on both
homophilic and heterophilic graphs. However, GPRGNN needs more supervised
information in order to learn optimal GPR weights. On the contrary, pGNNs need
less supervised information to obtain similar results because $\Theta^{(2)}$
acts like a hyperplane for classification. pGNNs could work better under weak
supervised information. Our analysis is also verified by the experimental
results in Sec. 5.

We also provide an upper-bounding risk of pGNNs by Theorem 4 in Sec. C.1 to
study the effect of the hyperparameter $\mu$ on the performance of pGNNs.
Theorem 4 shows that the risk of pGNNs is upper-bounded by the sum of three
terms: the risk of label prediction using only the original node features
${\mathbf{X}}$, the norm of $p$-Laplacian diffusion on ${\mathbf{X}}$, and the
magnitude of the noise in ${\mathbf{X}}$. $\mu$ controls the trade-off between
these three terms. The smaller $\mu$, the more weights on the $p$-Laplacian
diffusion term and the noise term and the less weights on the the other term
and vice versa.

##  4 Spectral Views of $p$-Laplacian Message Passing

In this section, we theoretically demonstrate that $p$-Laplacian message
passing is an approximation of a polynomial graph filter defined on the
spectral domain of $p$-Laplacian. We show by spectral analysis that
$p$-Laplacian message passing works as low-high-pass filters.

$p$-Eigenvalues and $p$-Eigenvectors of the Graph $p$-Laplacian. We first
introduce the definitions of $p$-eigenvalues and $p$-eigenvectors of
$p$-Laplacian. Let $\phi_{p}:\mathbb{R}\rightarrow\mathbb{R}$ defined as
$\phi_{p}(u)=\|u\|^{p-2}u,\text{ for }u\in\mathbb{R},u\neq 0$. Note that
$\phi_{2}(u)=u$. For notational simplicity, we denote by
$\phi_{p}({\mathbf{u}})=(\phi_{p}({u}_{1}),\dots,\phi_{p}({u}_{N}))^{\top}$
for ${\mathbf{u}}\in\mathbb{R}^{N}$ and
$\Phi_{p}({\mathbf{U}})=(\phi_{p}({\mathbf{U}}_{:,1}),\dots,\phi_{p}({\mathbf{U}}_{:,N}))$
for ${\mathbf{U}}\in\mathbb{R}^{N\times N}$ and
${\mathbf{U}}_{:,i}\in\mathbb{R}^{N}$ is the $i^{th}$ column vector of
${\mathbf{U}}$.

######  Definition 4 ($p$-Eigenvector and $p$-Eigenvalue).

A vector ${\mathbf{u}}\in\mathbb{R}^{N}$ is a $p$-eigenvector of $\Delta_{p}$
if it satisfies the equation

$\left(\Delta_{p}{\mathbf{u}}\right)_{i}=\lambda\phi_{p}({u}_{i}),\quad\text{for
all }i\in[N],$

where $\lambda\in\mathbb{R}$ is a real value referred as a $p$-eigenvalue of
$\Delta_{p}$ associated with the $p$-eigenvector ${\mathbf{u}}$.

######  Definition 5 ($p$-Orthogonal (Luo, Huang, Ding, and Nie, 2010)).

Given two vectors ${\mathbf{u}},{\mathbf{v}}\in\mathbb{R}^{N}$ with
${\mathbf{u}},{\mathbf{v}}\neq\bm{0}$, we call that ${\mathbf{u}}$ and
${\mathbf{v}}$ is $p$-orthogonal if

$\phi_{p}({\mathbf{u}})^{\top}\phi_{p}({\mathbf{v}})=\sum_{i=1}^{N}\phi_{p}({u}_{i})\phi_{p}(v_{i})=0.$

Luo et al. (2010) demonstrated that the $p$-eigenvectors of $\Delta_{p}$ are
$p$-orthogonal to each other (see Theorem 5 in Sec. C.2 for details).
Therefore, the space spanned by the multiple $p$-eigenvectors of $\Delta_{p}$
is $p$-orthogonal. Additionally, we demonstrate that the $p$-eigen-
decomposition of $\Delta_{p}$ is given by:
$\Delta_{p}=\Phi_{p}({\mathbf{U}}){\mathbf{\Lambda}}\Phi_{p}({\mathbf{U}})^{\top}$
(see Theorem 6 in Sec. C.3 for details), where ${\mathbf{U}}$ is a matrix of
$p$-eigenvectors of $\Delta_{p}$ and ${\mathbf{\Lambda}}$ is a diagonal matrix
in which the diagonal is the $p$-eigenvalues of $\Delta_{p}$.

Graph Convolutions based on $p$-Laplacian. Based on Theorem 5, the graph
Fourier Transform $\hat{f}$ of any function $f$ on the vertices of
${\mathcal{G}}$ can be defined as the expansion of $f$ in terms of
$\Phi({\mathbf{U}})$ where ${\mathbf{U}}$ is the matrix of $p$-eigenvectors of
$\Delta_{p}$: $\hat{f}=\Phi_{p}({\mathbf{U}})^{\top}f$. Similarly, the inverse
graph Fourier transform is then given by: $f=\Phi_{p}({\mathbf{U}})\hat{f}$.
Therefore, a signal ${\mathbf{X}}\in\mathbb{R}^{N\times c}$ being filtered by
a spectral filter $g_{\theta}$ can be expressed formally as:
$g_{\bm{\theta}}\star{\mathbf{X}}=\Phi_{p}({\mathbf{U}})\hat{g}_{\bm{\theta}}({\mathbf{\Lambda}})\Phi_{p}({\mathbf{U}})^{\top}{\mathbf{X}}$,
where ${\mathbf{\Lambda}}$ denotes a diagonal matrix in which the diagonal
corresponds to the $p$-eigenvalues $\\{\lambda_{l}\\}_{l=0,\dots,N-1}$ of
$\Delta_{p}$ and $\hat{g}_{\bm{\theta}}({\mathbf{\Lambda}})$ denotes a
diagonal matrix in which the diagonal corresponds to spectral filter
coefficients. Let $\hat{g}_{\bm{\theta}}$ be a polynomial filter defined as
$\hat{g}_{\bm{\theta}}=\sum_{k=0}^{K-1}\theta_{k}\lambda_{l}^{k}$, where the
parameter
$\bm{\theta}=[\theta_{0},\dots,\theta_{K-1}]^{\top}\in\mathbb{R}^{K}$ is a
vector of polynomial coefficients. By the $p$-eigen-decomposition of
$p$-Laplacian, we have

$g_{\bm{\theta}}\star{\mathbf{X}}\approx\sum_{k=0}^{K-1}\theta_{k}\Phi_{p}({\mathbf{U}}){\mathbf{\Lambda}}^{k}\Phi_{p}({\mathbf{U}})^{\top}{\mathbf{X}}=\sum_{k=0}^{K-1}\theta_{k}\Delta_{p}^{k}{\mathbf{X}}$
(18)

######  Theorem 3.

The $K$-step $p$-Laplacian message passing is a $K$-order polynomial
approximation to the graph filter given by Eq. 18.

Proof see Sec. D.3. Theorem 3 indicates that $p$-Laplacian message passing is
implicitly a polynomial spectral filter defined on the spectral domain of
$p$-Laplacian.

Spectral Analysis of $p$-Laplacian Message Passing. Here, we analyze the
spectral propecties of $p$-Laplacian message passing. We can approximately
view $p$-Laplacian message pasing as a filter of a linear combination of $K$
spectral filters
$g({\mathbf{\Lambda}})^{(0)},g({\mathbf{\Lambda}})^{(1)},\dots,g({\mathbf{\Lambda}})^{(K-1)}$
with each spectral filter defined to be
$g({\mathbf{\Lambda}})^{(k)}:=(\bm{\alpha}{\mathbf{D}}^{-1/2}{\mathbf{M}}{\mathbf{D}}^{-1/2})^{k}$
where
${M}_{i,j}={W}_{i,j}\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}\|^{p-2}$
for $i,j\in[N]$ and ${\mathbf{F}}$ is the matrix of node embeddings. We can
study the properties of $p$-Laplacian message passing by studying the spectral
properties of $\bm{\alpha}{\mathbf{D}}^{-1/2}{\mathbf{M}}{\mathbf{D}}^{-1/2}$
as given below.

######  Proposition 1.

Given a connected graph
${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$ with node
embeddings ${\mathbf{F}}$ and the $p$-Laplacian $\Delta_{p}$ with its
$p$-eigenvectors $\\{{\mathbf{u}}^{(l)}\\}_{l=0,1,\dots,N-1}$ and the
$p$-eigenvalues $\\{\lambda_{l}\\}_{l=0,1,\dots,N-1}$. Let
$g_{p}(\lambda_{i-1}):=\alpha_{i,i}\sum_{j}{D}_{i,i}^{-1/2}{M}_{i,j}{D}_{j,j}^{-1/2}$
for $i\in[N]$ be the filters defined on the spectral domain of $\Delta_{p}$,
where ${M}_{i,j}={W}_{i,j}\|\nabla f([i,j])\|^{p-2}$, $(\nabla f)([i,j])$ is
the graph gradient of the edge between node $i$ and $j$ and $\|\nabla f(i)\|$
is the norm of graph gradient at $i$. $N_{i}$ denotes the number of edges
connected to $i$, $N_{\text{min}}=\min\\{N_{j}\\}_{j\in[N]}$, and
$k=\operatorname*{arg\,max}_{j}(\\{|{u}_{j}^{(l)}|/\sqrt{{D}_{l,l}}\\}_{j\in[N];l=0,\dots,N-1})$,
then

  1. 1.

When $p=2$, $g_{p}(\lambda_{i-1})$ works as low-high-pass filters.

  2. 2.

When $p>2$, if $\|\nabla f(i)\|\leq 2^{(p-1)/(p-2)}$, $g_{p}(\lambda_{i-1})$
works as low-high-pass filters on node $i$ and $g_{p}(\lambda_{i-1})$ works as
low-pass filters on $i$ when $\|\nabla f(i)\|\geq 2^{(p-1)/(p-2)}$.

  3. 3.

When $1\leq p<2$, if $0\leq\|\nabla f(i)\|\leq 2(2\sqrt{N_{k}})^{1/(p-2)}$,
$g_{p}(\lambda_{i-1})$ works as low-pass filters on node $i$ and
$g_{p}(\lambda_{i-1})$ works as low-high-pass filters on $i$ when $\|\nabla
f(i)\|\geq 2\left(2\sqrt{N_{k}}\right)^{1/(p-2)}$. Specifically, when $p=1$,
$N_{k}$ can be replaced by $N_{\text{min}}$.

Proof see Sec. D.7. Proposition 1 shows that when $p\neq 2$, $p$-Laplacian
message passing adaptively works as low-pass or low-high-pass filters on node
$i$ in terms of the degree of local node embedding variation around $i$, i.e.
the norm of the graph gradient $\|\nabla f(i)\|$ at node $i$. When $p=2$,
$p$-Laplacian message passing works as low-high-pass filters on node $i$
regardless of the value of $\|\nabla f(i)\|$. When $p>2$, $p$-Laplacian
message passing works as low-pass filters on node $i$ for large $\|\nabla
f(i)\|$ and works as low-high-pass filters for small $\|\nabla f(i)\|$.
Therefore, pGNNs with $p>2$ can work very effectively on graphs with strong
homophily. When $1\leq p<2$, $p$-Laplacian message passing works as low-pass
filters for small $\|\nabla f(i)\|$ and works as low-high-pass filters for
large $\|\nabla f(i)\|$. Thus, pGNNs with $1\leq p<2$ can work effectively on
graphs with low homophily, i.e. heterophilic graphs.

######  Remark 4.

We say that a filter works as a low-pass filter if the low frequencies
dominate the other frequencies and it works as a low-high-pass filter if the
low frequencies and high frequencies dominate the middle frequencies, i.e.,
the frequencies around the cutoff frequency.

######  Remark 5.

Previous works (Wu et al., 2019; Klicpera et al., 2019) have shown that GCN,
SGC, APPNP work as low-pass filters. The reason is that they have added the
self-loop to the adjacency matrix, which will shrink the spectral range of the
Laplacian from $[0,2]$ to approximately $[0,1.5]$ and causes them work as low-
pass filters (Wu et al., 2019). On the contrary, we did not add self-loop to
the adjacancy matrix and therefore the spectral filter work as low-high-pass
filters in our case for $p=2$.

##  5 Empirical Studies

Here, we empirically study the effectiveness of pGNNs for semi-supervised node
classification using and real-world benchmark and synthetic datasets with
heterophily and strong homophily. The experimental results are also used to
validate our theoretical findings presented previously.

Table 1: Heterophilious results. Averaged accuracy (%) for 100 runs. Best
results outlined in bold and the results within $95\%$ confidence interval of
the best results are outlined in underlined bold.

Method | Chameleon | Squirrel | Actor | Wisconsin | Texas | Cornell
---|---|---|---|---|---|---
MLP | $48.02_{\pm 1.72}$ | $\bm{33.80}_{\pm 1.05}$ | $39.68_{\pm 1.43}$ | $93.56_{\pm 3.14}$ | $79.50_{\pm 10.62}$ | $80.30_{\pm 11.38}$
GCN | $34.54_{\pm 2.78}$ | $25.28_{\pm 1.55}$ | $31.28_{\pm 2.04}$ | $61.93_{\pm 3.00}$ | $56.54_{\pm 17.02}$ | $51.36_{\pm 4.59}$
SGC | $34.76_{\pm 4.55}$ | $25.49_{\pm 1.63}$ | $30.98_{\pm 3.80}$ | $66.94_{\pm 2.58}$ | $59.99_{\pm 9.95}$ | $44.39_{\pm 5.88}$
GAT | $45.16_{\pm 2.10}$ | $31.41_{\pm 0.98}$ | $34.11_{\pm 1.28}$ | $65.64_{\pm 6.29}$ | $56.41_{\pm 13.01}$ | $43.94_{\pm 7.33}$
JKNet | $33.28_{\pm 3.59}$ | $25.82_{\pm 1.58}$ | $29.77_{\pm 2.61}$ | $61.08_{\pm 3.71}$ | $59.65_{\pm 12.62}$ | $55.34_{\pm 4.43}$
APPNP | $36.18_{\pm 2.81}$ | $26.85_{\pm 1.48}$ | $31.26_{\pm 2.52}$ | $64.59_{\pm 3.49}$ | $82.90_{\pm 5.08}$ | $66.47_{\pm 9.34}$
GPRGNN | $43.67_{\pm 2.27}$ | $31.27_{\pm 1.76}$ | $36.63_{\pm 1.22}$ | $88.54_{\pm 4.94}$ | $80.74_{\pm 6.76}$ | $78.95_{\pm 8.52}$
1.0GNN | $\bm{48.86}_{\pm 1.95}$ | $\bm{\underline{33.75}}_{\pm 1.50}$ | $\bm{40.62}_{\pm 1.25}$ | $\bm{95.37}_{\pm 2.06}$ | $84.06_{\pm 7.41}$ | $\bm{82.16}_{\pm 8.62}$
1.5GNN | $\bm{\underline{48.74}}_{\pm 1.62}$ | $33.33_{\pm 1.45}$ | $\bm{\underline{40.35}}_{\pm 1.35}$ | $\bm{\underline{95.24}}_{\pm 2.01}$ | $84.46_{\pm 7.79}$ | $78.47_{\pm 6.87}$
2.0GNN | $\bm{\underline{48.77}}_{\pm 1.87}$ | $33.60_{\pm 1.47}$ | $40.07_{\pm 1.17}$ | $91.15_{\pm 2.76}$ | $\bm{87.96}_{\pm 6.27}$ | $72.04_{\pm 8.22}$
2.5GNN | $\bm{\underline{48.80}}_{\pm 1.77}$ | $\bm{\underline{33.79}}_{\pm 1.45}$ | $39.80_{\pm 1.31}$ | $87.08_{\pm 2.69}$ | $83.01_{\pm 6.80}$ | $70.31_{\pm 8.84}$

Datasets and Experimental Setup. We use seven homophilic benchmark datasets:
citation graphs Cora, CiteSeer, PubMed (Sen et al., 2008), Amazon co-purchase
graphs Computers, Photo, coauthor graphs CS, Physics (Shchur et al., 2018),
and six heterophilic benchmark datasets: Wikipedia graphs Chameleon, Squirrel
(Rozemberczki et al., 2021), the Actor co-occurrence graph, webpage graphs
Wisconsin, Texas, Cornell (Pei et al., 2020). The node classification tasks
are conducted in the transductive setting. Following Chien et al. (2021), we
use the sparse splitting ($2.5\%/2.5\%/95\%$) and the dense splitting
($60\%/20\%/20\%$) to randomly split the homophilic and heterophilic graphs
into training$/$validation$/$testing sets, respectively. Note that we directly
used the data from Pytorch Geometric library (Fey & Lenssen, 2019) where they
did not transform Chameleon and Squirrel to undirected graphs, which is
different from Chien et al. (2021) where they did so. Dataset statistics and
their levels of homophily are presented in Appendix E.

Baselines. We compare pGNN with seven models, including MLP, GCN (Kipf &
Welling, 2017), SGC (Wu et al., 2019), GAT (Velickovic et al., 2018), JKNet
(Xu et al., 2018), APPNP (Klicpera et al., 2019), GPRGNN (Chien et al., 2021).
We use the Pytorch Geometric library to implement all baselines except GPRGNN.
For GPRGNN, we use the code released by the
authors444https://github.com/jianhao2016/GPRGNN. The details of hyperparameter
settings are deferred to Sec. E.3.

Figure 1: Aggregation weight entropy distribution of graphs. Low entropy means
high degree of concentration, vice versa. An entropy of zero means all
aggregation weights are on one source node.

Superior Performance on Real-World Heterophilic Datasets. The results on
_homophilic_ benchmark datasets are deferred to Sec. F.1, which show that
pGNNs obtains competitive performance against state-of-the-art GNNs on
homophilic datasets. Table 1 summarizes the results on _heterophilic_
benchmark datasets. Table 1 shows that pGNNs significantly dominate the
baselines and 1.0GNN obtains the best performance on all heterophilic graphs
except the Texas dataset. For Texas, 2.0GNN is the best. We also observe that
MLP works very well and significantly outperforms most GNN baselines, which
indicates that the graph topology is not informative for label prediction on
these heterophilic graphs. Therefore, propagating and transforming node
features over the graph topology could lead to worse performance than MLP.
Unlike ordinary GNNs, pGNNs can adaptively learn aggregation weights and
ignore edges that are not informative for label prediction and thus could work
better. It confirms our theoretical findings presented in previous sections.
Note that GAT can also learn aggregation weights, i.e. the attention weights.
However, the aggregation weights learned by GAT are significantly distinct
from that of pGNNs, as we will show following.

Interpretability of the Learned Aggregation Weights of pGNNs. We showcase the
interpretability of the learned aggregation weights
$\alpha_{i,i}{D}_{i,i}^{-1/2}{M}_{i,j}{D}_{j,j}^{-1/2}$ of pGNNs by studying
its entropy distribution, along with the attention weights of GAT on real-
world datasets. Denote $\\{{A}_{i,j}\\}_{j\in{\mathcal{N}}_{i}}$ as the
aggregation weights of node $i$ and its neighbors. For GAT,
$\\{{A}_{i,j}\\}_{j\in{\mathcal{N}}_{i}}$ are referred as the attention
weights (in the first layer) and for pGNNs are
$\alpha_{i,i}{D}_{i,i}^{-1/2}{M}_{i,j}{D}_{j,j}^{-1/2}$. For any node $i$,
$\\{{A}_{i,j}\\}_{j\in{\mathcal{N}}_{i}}$ forms a discrete probability
distribution over all its neighbors with the entropy given by
$H(\\{{A}_{i,j}\\}_{j\in{\mathcal{N}}_{i}})=-\sum_{j\in{\mathcal{N}}_{i}}{A}_{i,j}\log({A}_{i,j})$.
Low entropy means high degree of concentration and vice versa. An entropy of
zero means all aggregation weights or attentions are on one source node. The
uniform distribution has the highest entropy of $log({D}_{i,i})$. Fig. 1
reports the results on Computers, Wisconsin and we defer more results on other
datasets to Sec. F.2 due to space limit. Fig. 1 shows that the aggregation
weight entropy distributions of GAT and pGNNs on Computers (homophily) are
both similar to the uniform case. It indicates the original graph topology of
Computers is very helpful for label prediction and therefore GNNs could work
very well on Computers. However, for Wisconsin (heterophily), the entropy
distribution of pGNNs is significantly different from that of GAT and the
uniform case. Most entropy of pGNNs is around zero, which means that most
aggregation weights are on one source node. It states that the original graph
topology of Wisconsin is not helpful for label prediction, which explains why
MLP works well on Wisconsin. On the contrary, the entropy distribution of GAT
is similar to the uniform case and therefore GAT works similarly to GCN and is
significantly worse than pGNNs on Wisconsin. Similar results can be observed
on the experiments on more datasets in Sec. F.2.

Figure 2: Averaged accuracy (%) on cSBM (sparse split) for 20 runs. Best view
in colors.

Results on cSBM Datasets. We examine the performance of pGNNs on heterophilic
graphs whose topology is informative for label prediction using synthetic
graphs generated by cSBM (Deshpande et al., 2018) with
$\phi\in\\{-1,-0.75,\dots,1\\}$. We use the same settings of cSBM used in
Chien et al. (2021). Due to the space limit, we refer the readers to Chien et
al. (2021) for more details of cSBM dataset. Fig. 2 reports the results on
cSBM using sparse splitting (for results on cSBM with dense splitting see Sec.
F.3). Fig. 2 shows that when $\phi\leq-0.5$ (heterophilic graphs), 2.0GNN
obtains the best performance and pGNNs and GPRGNN significantly dominate the
others. It validates the effectiveness of pGNNs on heterophilic graphs.
Moreover, 2.0GNN works better than GPRGNN and it again confirms that 2.0GNN is
more superior under weak supervision (2.5% training rate), as stated in Remark
3. Note that 1.0GNN and 1.5GNN are not better than 2.0GNN, the reason could be
the iteration algorithms Eq. 11 with $p=1,1.5$ are not as stable as the one
with $p=2$. When the graph topology is almost non-informative for label
prediction ($\phi=-0.25,0$), The performance of pGNNs is close to MLP and they
outperform the other baselines. Again, it validates that pGNNs can erase non-
informative edges and work as well as MLP and confirms the statements in
Theorem 4. When the graph is homophilic ($\phi\geq 0.25$), 1.5GNN is the best
on weak homophilic graphs ($\phi=0.25,0.5$) and pGNNs work competitively with
all GNN baselines on strong homophilic graphs ($\phi\geq 0.75$).

Figure 3: Averaged accuracy (%) on graphs with noisy edges for 20 runs. Best
view in colors.

Results on Datasets with Noisy Edges. We conduct experiments to evaluate the
performance of pGNNs on graphs with noisy edges by randomly adding edges to
the graphs and randomly remove the same number of original edges. We define
the random edge rate as $r:=\frac{\\#\text{random edges}}{\\#\text{all
edges}}$. The experiments are conducted on 4 homophilic datasets (Computers,
Photo, CS, Physics) and 2 heterophilic datasets (Wisconsin, Texas) with
$r=0.25,0.5,1$. Fig. 3 reports the results on Computers, Wisconsin and we
defer more results to Sec. F.4. Fig. 3 shows that pGNNs significantly
outperform all baselines. Specifically, 1.5GNN obtains the best performance on
Computers, and 1.5GNN and 2.0GNN even work as well as MLP on Computers with
completely random edges ($r=1$). For Wisconsin, 1.0GNN is the best, and 1.0GNN
and 1.5GNN significantly dominate the others. We also observed that APPNP and
GPRGNN, whose architectures are analogical to 2.0GNN, also work better than
other GNNs. Nevertheless, they are significantly outperformed by pGNNs
overall. Similar results can be observed in the experiments conducted on more
datasets as presented in Sec. F.4.

##  6 Conclusion

We have addressed the problem of generalizing GNNs to heterophilic graphs and
graphs with noisy edges. To this end, we derived a novel $p$-Laplacian message
passing scheme from a discrete regularization framework and proposed a new
pGNN architecture. We theoretically demonstrate our method works as low-high-
pass filters and thereby applicable to both homophilic and heterophilic
graphs. We empirically validate our theoretical results and show the
advantages of our methods on heterophilic graphs and graphs with non-
informative topologies.

Like most existing spectral based GNN models, e.g., GCN (Kipf & Welling,
2017), SGC (Wu et al., 2019), the main restriction of pGNNs is the relatively
high space cost compared to GraphSage (Hamilton et al., 2017), especially for
extremely large graphs. Integrating pGNNs and $p$-Laplacian message passing
with some popular subgraph sampling techniques so that pGNNs (or its variants)
could scale to large graphs would be an interesting future work. We refer the
reader to Appendix B for further discussions on the potential extensions of
pGNNs.

## Acknowledgements

Guoji would like to thank Prof. Xin Yao and Dr. Yunwen Lei for their sincere
and selfless helps and supports.

## Reproducibility Statement

In order to ensure reproducibility, we have made the efforts in the following
respects: (1) Provide our code as the supplementary material; (2) Provide
self-contained proofs of the main claims in Appendices C and D; (3) Provide
more details on experimental configurations in Appendix E and experimental
results in Appendix F. All the datasets are publicly available as described in
the main text. The code of pGNNs is available at https://github.com/guoji-
fu/pGNNs.

## References

  * Abu-El-Haija et al. (2018) Abu-El-Haija, S., Perozzi, B., Al-Rfou, R., and Alemi, A. A.  Watch your step: Learning node embeddings via graph attention.  In Bengio, S., Wallach, H. M., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 31, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pp. 9198–9208, 2018.
  * Abu-El-Haija et al. (2019) Abu-El-Haija, S., Kapoor, A., Perozzi, B., and Lee, J.  N-GCN: multi-scale graph convolution for semi-supervised node classification.  In Globerson, A. and Silva, R. (eds.), _Proceedings of the Thirty-Fifth Conference on Uncertainty in Artificial Intelligence, UAI 2019, Tel Aviv, Israel, July 22-25, 2019_ , volume 115 of _Proceedings of Machine Learning Research_ , pp. 841–851. AUAI Press, 2019.
  * Amghibech (2003) Amghibech, S.  Eigenvalues of the discrete p-laplacian for graphs.  _Ars Combinatoria_ , 67:283–302, 2003.
  * Atwood & Towsley (2016) Atwood, J. and Towsley, D.  Diffusion-convolutional neural networks.  In Lee, D. D., Sugiyama, M., von Luxburg, U., Guyon, I., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 29, NeurIPS 2016, December 5-10, 2016, Barcelona, Spain_ , pp. 1993–2001, 2016\.
  * Bach & Jordan (2004) Bach, F. and Jordan, M.  Learning spectral clustering.  _Advances in Neural Information Processing Systems, NIPS 2004_ , 16(2):305–312, 2004.
  * Battaglia et al. (2018) Battaglia, P. W., Hamrick, J. B., Bapst, V., Sanchez-Gonzalez, A., Zambaldi, V. F., Malinowski, M., Tacchetti, A., Raposo, D., Santoro, A., Faulkner, R., Gülçehre, Ç., Song, H. F., Ballard, A. J., Gilmer, J., Dahl, G. E., Vaswani, A., Allen, K. R., Nash, C., Langston, V., Dyer, C., Heess, N., Wierstra, D., Kohli, P., Botvinick, M., Vinyals, O., Li, Y., and Pascanu, R.  Relational inductive biases, deep learning, and graph networks.  _CoRR_ , abs/1806.01261, 2018.
  * Belkin & Niyogi (2008) Belkin, M. and Niyogi, P.  Towards a theoretical foundation for laplacian-based manifold methods.  _Journal of Computer and System Sciences_ , 74(8):1289–1308, 2008.
  * Belkin et al. (2004) Belkin, M., Matveeva, I., and Niyogi, P.  Regularization and semi-supervised learning on large graphs.  In _The 17th Annual Conference on Learning Theory, COLT 2004, Banff, Canada, July 1-4, 2004_ , volume 3120, pp. 624–638, 2004.
  * Belkin et al. (2005) Belkin, M., Niyogi, P., and Sindhwani, V.  On manifold regularization.  In Cowell, R. G. and Ghahramani, Z. (eds.), _Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, AISTATS 2005, Bridgetown, Barbados, January 6-8, 2005_. Society for Artificial Intelligence and Statistics, 2005.
  * Bruna et al. (2014) Bruna, J., Zaremba, W., Szlam, A., and LeCun, Y.  Spectral networks and locally connected networks on graphs.  In _The 2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014_ , 2014.
  * Bühler & Hein (2009) Bühler, T. and Hein, M.  Spectral clustering based on the graph _p_ -laplacian.  In Danyluk, A. P., Bottou, L., and Littman, M. L. (eds.), _Proceedings of the 26th Annual International Conference on Machine Learning, ICML 2009, Montreal, Quebec, Canada, June 14-18, 2009_ , volume 382 of _ACM International Conference Proceeding Series_ , pp. 81–88. ACM, 2009.
  * Chen et al. (2018) Chen, J., Ma, T., and Xiao, C.  Fastgcn: Fast learning with graph convolutional networks via importance sampling.  In _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings_ , 2018.
  * Chien et al. (2021) Chien, E., Peng, J., Li, P., and Milenkovic, O.  Adaptive universal generalized pagerank graph neural network.  In _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net, 2021.
  * Defferrard et al. (2016) Defferrard, M., Bresson, X., and Vandergheynst, P.  Convolutional neural networks on graphs with fast localized spectral filtering.  In _Advances in Neural Information Processing Systems 29, NeurIPS 2016, December 5-10, 2016, Barcelona, Spain_ , pp. 3837–3845, 2016\.
  * Deshpande et al. (2018) Deshpande, Y., Sen, S., Montanari, A., and Mossel, E.  Contextual stochastic block models.  In Bengio, S., Wallach, H. M., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 31, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pp. 8590–8602, 2018.
  * Duvenaud et al. (2015) Duvenaud, D., Maclaurin, D., Aguilera-Iparraguirre, J., Gómez-Bombarelli, R., Hirzel, T., Aspuru-Guzik, A., and Adams, R. P.  Convolutional networks on graphs for learning molecular fingerprints.  In Cortes, C., Lawrence, N. D., Lee, D. D., Sugiyama, M., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 28, NeurIPS 2015 December 7-12, 2015, Montreal, Quebec, Canada_ , pp. 2224–2232, 2015.
  * Fey & Lenssen (2019) Fey, M. and Lenssen, J. E.  Fast graph representation learning with pytorch geometric.  _CoRR_ , abs/1903.02428, 2019.
  * Fu et al. (2020) Fu, G., Hou, Y., Zhang, J., Ma, K., Kamhoua, B. F., and Cheng, J.  Understanding graph neural networks from graph signal denoising perspectives.  _CoRR_ , abs/2006.04386, 2020.
  * Gama et al. (2020) Gama, F., Bruna, J., and Ribeiro, A.  Stability properties of graph neural networks.  _IEEE Transactions on Signal Processing_ , 68:5680–5695, 2020.
  * Garg et al. (2020) Garg, V. K., Jegelka, S., and Jaakkola, T. S.  Generalization and representational limits of graph neural networks.  In _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 3419–3430. PMLR, 2020\.
  * Grandvalet & Bengio (2004) Grandvalet, Y. and Bengio, Y.  Semi-supervised learning by entropy minimization.  In _Advances in Neural Information Processing Systems 17, NIPS 2004, December 13-18, 2004, Vancouver, British Columbia, Canada]_ , pp. 529–536, 2004.
  * Hamilton et al. (2017) Hamilton, W. L., Ying, Z., and Leskovec, J.  Inductive representation learning on large graphs.  In _Advances in Neural Information Processing Systems 30, NeurIPS 2017, 4-9 December 2017, Long Beach, CA, USA_ , pp. 1024–1034, 2017\.
  * Hein (2006) Hein, M.  Uniform convergence of adaptive graph-based regularization.  In _The 19th Annual Conference on Learning Theory, COLT 2006, Pittsburgh, PA, USA, June 22-25, 2006_ , volume 4005, pp. 50–64, 2006.
  * Henaff et al. (2015) Henaff, M., Bruna, J., and LeCun, Y.  Deep convolutional networks on graph-structured data.  _CoRR_ , abs/1506.05163, 2015.
  * Hu et al. (2020) Hu, W., Fey, M., Zitnik, M., Dong, Y., Ren, H., Liu, B., Catasta, M., and Leskovec, J.  Open graph benchmark: Datasets for machine learning on graphs.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M., and Lin, H. (eds.), _Advances in Neural Information Processing Systems 33, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Huang et al. (2021) Huang, Q., He, H., Singh, A., Lim, S., and Benson, A. R.  Combining label propagation and simple models out-performs graph neural networks.  In _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net, 2021.
  * Kipf & Welling (2017) Kipf, T. N. and Welling, M.  Semi-supervised classification with graph convolutional networks.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_ , 2017.
  * Klicpera et al. (2019) Klicpera, J., Bojchevski, A., and Günnemann, S.  Predict then propagate: Graph neural networks meet personalized pagerank.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Li et al. (2019) Li, G., Müller, M., Thabet, A. K., and Ghanem, B.  Deepgcns: Can gcns go as deep as cnns?  In _2019 IEEE/CVF International Conference on Computer Vision, ICCV 2019, Seoul, Korea (South), October 27 - November 2, 2019_ , pp. 9266–9275. IEEE, 2019.
  * Li et al. (2018) Li, Q., Han, Z., and Wu, X.  Deeper insights into graph convolutional networks for semi-supervised learning.  In _Thirty-Second AAAI Conference on Artificial Intelligence, AAAI 2018, New Orleans, Louisiana, USA, February 2-7, 2018_ , pp. 3538–3545, 2018.
  * Liao et al. (2019) Liao, R., Zhao, Z., Urtasun, R., and Zemel, R. S.  Lanczosnet: Multi-scale deep graph convolutional networks.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Liu et al. (2014) Liu, F., Chakraborty, S., Li, F., Liu, Y., and Lozano, A. C.  Bayesian regularization via graph laplacian.  _Bayesian Analysis_ , 9(2):449–474, 2014.
  * Liu & Han (2013) Liu, J. and Han, J.  Spectral clustering.  In Aggarwal, C. C. and Reddy, C. K. (eds.), _Data Clustering: Algorithms and Applications_ , pp. 177–200. CRC Press, 2013.
  * Liu et al. (2021) Liu, X., Jin, W., Ma, Y., Li, Y., Liu, H., Wang, Y., Yan, M., and Tang, J.  Elastic graph neural networks.  In Meila, M. and Zhang, T. (eds.), _Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event_ , volume 139 of _Proceedings of Machine Learning Research_ , pp. 6837–6849. PMLR, 2021.
  * Loukas (2020) Loukas, A.  What graph neural networks cannot learn: depth vs width.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Luo et al. (2010) Luo, D., Huang, H., Ding, C. H. Q., and Nie, F.  On the eigenvectors of _p_ -laplacian.  _Machine Learning_ , 81(1):37–51, 2010.
  * Marcheggiani & Titov (2017) Marcheggiani, D. and Titov, I.  Encoding sentences with graph convolutional networks for semantic role labeling.  In Palmer, M., Hwa, R., and Riedel, S. (eds.), _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing, EMNLP 2017, Copenhagen, Denmark, September 9-11, 2017_ , pp. 1506–1515. Association for Computational Linguistics, 2017.
  * Nadler et al. (2009) Nadler, B., Srebro, N., and Zhou, X.  Semi-supervised learning with the graph laplacian: The limit of infinite unlabelled data.  _Advances in Neural Information Processing Systems, NIPS 2009_ , 22:1330–1338, 2009.
  * Niepert et al. (2016) Niepert, M., Ahmed, M., and Kutzkov, K.  Learning convolutional neural networks for graphs.  In Balcan, M. and Weinberger, K. Q. (eds.), _Proceedings of the 33nd International Conference on Machine Learning, ICML 2016, New York City, NY, USA, June 19-24, 2016_ , volume 48 of _JMLR Workshop and Conference Proceedings_ , pp. 2014–2023. JMLR.org, 2016.
  * Niyogi (2013) Niyogi, P.  Manifold regularization and semi-supervised learning: some theoretical analyses.  _Journal of Machine Learning Research_ , 14(1):1229–1250, 2013.
  * NT & Maehara (2019) NT, H. and Maehara, T.  Revisiting graph neural networks: All we have is low-pass filters.  _CoRR_ , abs/1905.09550, 2019.
  * Oono & Suzuki (2020) Oono, K. and Suzuki, T.  Graph neural networks exponentially lose expressive power for node classification.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Page et al. (1999) Page, L., Brin, S., Motwani, R., and Winograd, T.  The pagerank citation ranking: Bringing order to the web.  Technical report, Stanford InfoLab, 1999.
  * Pei et al. (2020) Pei, H., Wei, B., Chang, K. C., Lei, Y., and Yang, B.  Geom-gcn: Geometric graph convolutional networks.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Rozemberczki et al. (2021) Rozemberczki, B., Allen, C., and Sarkar, R.  Multi-scale attributed node embedding.  _Journal of Complex Networks_ , 9(2), 2021.
  * Satorras & Estrach (2018) Satorras, V. G. and Estrach, J. B.  Few-shot learning with graph neural networks.  In _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings_. OpenReview.net, 2018.
  * Sen et al. (2008) Sen, P., Namata, G., Bilgic, M., Getoor, L., Gallagher, B., and Eliassi-Rad, T.  Collective classification in network data.  _AI Magazine_ , 29(3):93–106, 2008.
  * Shchur et al. (2018) Shchur, O., Mumme, M., Bojchevski, A., and Günnemann, S.  Pitfalls of graph neural network evaluation.  _CoRR_ , abs/1811.05868, 2018.
  * Sindhwani et al. (2005) Sindhwani, V., Niyogi, P., Belkin, M., and Keerthi, S.  Linear manifold regularization for large scale semi-supervised learning.  In _Proceedings of the 22nd ICML Workshop on Learning with Partially Classified Training Data_ , volume 28, 2005.
  * Slepcev & Thorpe (2017) Slepcev, D. and Thorpe, M.  Analysis of $p$-laplacian regularization in semi-supervised learning.  _CoRR_ , abs/1707.06213, 2017.
  * Smola & Kondor (2003) Smola, A. J. and Kondor, R.  Kernels and regularization on graphs.  In Schölkopf, B. and Warmuth, M. K. (eds.), _Computational Learning Theory and Kernel Machines, 16th Annual Conference on Computational Learning Theory and 7th Kernel Workshop, COLT/Kernel 2003, Washington, DC, USA, August 24-27, 2003, Proceedings_ , volume 2777 of _Lecture Notes in Computer Science_ , pp. 144–158. Springer, 2003.
  * Thekumparampil et al. (2018) Thekumparampil, K. K., Wang, C., Oh, S., and Li, L.  Attention-based graph neural network for semi-supervised learning.  _CoRR_ , abs/1803.03735, 2018.
  * van der Maaten & Hinton (2008) van der Maaten, L. and Hinton, G.  Visualizing data using t-sne.  _Journal of Machine Learning Research_ , 9(86):2579–2605, 2008.
  * van Engelen & Hoos (2020) van Engelen, J. E. and Hoos, H. H.  A survey on semi-supervised learning.  _Machine Learning_ , 109(2):373–440, 2020.
  * Velickovic et al. (2018) Velickovic, P., Cucurull, G., Casanova, A., Romero, A., Liò, P., and Bengio, Y.  Graph attention networks.  In _6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018_ , 2018.
  * Velickovic et al. (2019) Velickovic, P., Fedus, W., Hamilton, W. L., Liò, P., Bengio, Y., and Hjelm, R. D.  Deep graph infomax.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Verma & Zhang (2019) Verma, S. and Zhang, Z.  Stability and generalization of graph convolutional neural networks.  In _The 25th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD 2019, Anchorage, AK, USA, August 4-8, 2019_ , pp. 1539–1548, 2019.
  * von Luxburg (2007) von Luxburg, U.  A tutorial on spectral clustering.  _Statistics and Computing_ , 17(4):395–416, 2007\.
  * Wang & Leskovec (2020) Wang, H. and Leskovec, J.  Unifying graph convolutional neural networks and label propagation.  _CoRR_ , abs/2002.06755, 2020.
  * Wu et al. (2019) Wu, F., Jr., A. H. S., Zhang, T., Fifty, C., Yu, T., and Weinberger, K. Q.  Simplifying graph convolutional networks.  In _Proceedings of the 36th International Conference on Machine Learning, ICML 2019, 9-15 June 2019, Long Beach, California, USA_ , pp. 6861–6871, 2019.
  * Wu et al. (2021) Wu, Z., Pan, S., Chen, F., Long, G., Zhang, C., and Yu, P. S.  A comprehensive survey on graph neural networks.  _IEEE Transactions on Neural Networks and Learning System_ , 32(1):4–24, 2021.
  * Xinyi & Chen (2019) Xinyi, Z. and Chen, L.  Capsule graph neural network.  In _7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019_. OpenReview.net, 2019.
  * Xu et al. (2018) Xu, K., Li, C., Tian, Y., Sonobe, T., Kawarabayashi, K., and Jegelka, S.  Representation learning on graphs with jumping knowledge networks.  In _Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholmsmässan, Stockholm, Sweden, July 10-15, 2018_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 5449–5458. PMLR, 2018.
  * Xu et al. (2019) Xu, K., Hu, W., Leskovec, J., and Jegelka, S.  How powerful are graph neural networks?  In _7th International Conference on Learning Representations, ICLR 2017 New Orleans, LA, USA, May 6-9, 2019, Conference Track Proceedings_ , 2019.
  * Ying et al. (2018) Ying, Z., You, J., Morris, C., Ren, X., Hamilton, W. L., and Leskovec, J.  Hierarchical graph representation learning with differentiable pooling.  In Bengio, S., Wallach, H. M., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R. (eds.), _Advances in Neural Information Processing Systems 31, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pp. 4805–4815, 2018.
  * Ying et al. (2019) Ying, Z., Bourgeois, D., You, J., Zitnik, M., and Leskovec, J.  Gnnexplainer: Generating explanations for graph neural networks.  In _Advances in Neural Information Processing Systems 32, NeurIPS 2019, 8-14 December 2019, Vancouver, BC, Canada_ , pp. 9240–9251, 2019\.
  * Zeng et al. (2020) Zeng, H., Zhou, H., Srivastava, A., Kannan, R., and Prasanna, V. K.  Graphsaint: Graph sampling based inductive learning method.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Zhou & Schölkopf (2005) Zhou, D. and Schölkopf, B.  Regularization on discrete spaces.  In _The 27th DAGM Symposium, Vienna, Austria, August 31 - September 2, 2005_ , volume 3663, pp. 361–368, 2005.
  * Zhou et al. (2003) Zhou, D., Bousquet, O., Lal, T. N., Weston, J., and Schölkopf, B.  Learning with local and global consistency.  In Thrun, S., Saul, L. K., and Schölkopf, B. (eds.), _Advances in Neural Information Processing Systems 16, NIPS 2003, December 8-13, 2003, Vancouver and Whistler, British Columbia, Canada]_ , pp. 321–328. MIT Press, 2003.
  * Zhou et al. (2020) Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C., and Sun, M.  Graph neural networks: A review of methods and applications.  _AI Open_ , 1:57–81, 2020.
  * Zhou & Belkin (2011) Zhou, X. and Belkin, M.  Semi-supervised learning by higher order regularization.  In Gordon, G. J., Dunson, D. B., and Dudík, M. (eds.), _Proceedings of the Fourteenth International Conference on Artificial Intelligence and Statistics, AISTATS 2011, Fort Lauderdale, USA, April 11-13, 2011_ , volume 15 of _JMLR Proceedings_ , pp. 892–900. JMLR.org, 2011.
  * Zhu et al. (2020) Zhu, J., Yan, Y., Zhao, L., Heimann, M., Akoglu, L., and Koutra, D.  Beyond homophily in graph neural networks: Current limitations and effective designs.  In _Advances in Neural Information Processing Systems 33, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Zhu et al. (2021) Zhu, J., Rossi, R. A., Rao, A., Mai, T., Lipka, N., Ahmed, N. K., and Koutra, D.  Graph neural networks with heterophily.  In _35th AAAI Conference on Artificial Intelligence, AAAI 2021, Virtual Event, February 2-9, 2021_ , pp. 11168–11176. AAAI Press, 2021\.
  * Zhu et al. (2003) Zhu, X., Ghahramani, Z., and Lafferty, J. D.  Semi-supervised learning using gaussian fields and harmonic functions.  In Fawcett, T. and Mishra, N. (eds.), _Machine Learning, Proceedings of the Twentieth International Conference, ICML 2003, August 21-24, 2003, Washington, DC, USA_ , pp. 912–919. AAAI Press, 2003.
  * Zitnik & Leskovec (2017) Zitnik, M. and Leskovec, J.  Predicting multicellular function through multi-layer tissue networks.  _Bioinform._ , 33(14):i190–i198, 2017.

Appendix

.tocmtappendix mtchapternone mtappendixsubsection

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries and Background
  3. 3 $p$-Laplacian based Graph Neural Networks
    1. 3.1 $p$-Laplacian Regularization Framework
    2. 3.2 $p$-Laplacian Message Passing and pGNN Architecture
  4. 4 Spectral Views of $p$-Laplacian Message Passing
  5. 5 Empirical Studies
  6. 6 Conclusion
  7. A Related Work
  8. B Discussions and Future Work
  9. C Additional Theorems
    1. C.1 Theorem 4 (Upper-Bounding Risk of pGNN)
    2. C.2 Theorem 5 ($p$-Orthogonal Theorem (Luo et al., 2010))
    3. C.3 Theorem 6 ($p$-Eigen-Decomposition of $\Delta_{p}$)
    4. C.4 Theorem 7 (Bounds of $p$-Eigenvalues)
  10. D Proof of Theorems
    1. D.1 Proof of Theorem 1
    2. D.2 Proof of Theorem 2
    3. D.3 Proof of Theorem 3
    4. D.4 Proof of Theorem 4
    5. D.5 Proof of Theorem 6
    6. D.6 Proof of Theorem 7
    7. D.7 Proof of Proposition 1
  11. E Dataset Statistics and Hyperparameters
    1. E.1 Illustration of Graph Gradient and Graph Divergence
    2. E.2 Dataset Statistics
    3. E.3 Hyperparameter Settings
  12. F Additional Experiments
    1. F.1 Experimental Results on Homophilic Benchmark Datasets
    2. F.2 Experimental Results of Aggregation Weight Entropy Distribution
    3. F.3 Experimental Results on cSBM
    4. F.4 Experimental Results on Graphs with Noisy Edges
    5. F.5 Experimental Results of Intergrating pGNNs with GCN and JKNet
    6. F.6 Experimental Results of pGNNs on PPI Dataset for Inductive Learning
    7. F.7 Experimental Results of pGNNs on OGBN arXiv Dataset
    8. F.8 Running Time of pGNNs
    9. F.9 Experimental Results on Benchmark Datasets for 64 Hidden Units
    10. F.10 Training Curves for $p=1$
    11. F.11 Visualization Results of Node Embeddings

##  Appendix A Related Work

Graph Neural Networks. Graph neural networks (GNNs) are a variant of neural
networks for graph-structured data, which can propagate and transform the node
features over the graph topology and exploit the information in the graphs.
Graph convolutional networks (GCNs) are one type of GNNs whose graph
convolution mechanisms or the message passing schemes were mainly inspired by
the field of graph signal processing. Bruna et al. (2014) defined a
nonparametric graph filter using the Fourier coefficients. Defferrard et al.
(2016) introduced Chebyshev polynomial to avoid computational expensive eigen-
decomposition of Laplacian and obtain localized spectral filters. GCN (Kipf &
Welling, 2017) used the first-order approximation and reparameterized trick to
simplify the spectral filters and obtain the layer-wise graph convolution. SGC
(Wu et al., 2019) further simplify GCN by removing non-linear transition
functions between each layer. Chen et al. (2018) propose importance sampling
to design an efficient variant of GCN. Xu et al. (2018) explored a jumping
knowledge architecture that flexibly leverages different neighborhood ranges
for each node to enable better structure-aware representation. Atwood &
Towsley (2016); Liao et al. (2019); Abu-El-Haija et al. (2019) exploited
multi-scale information by diffusing multi-hop neighbor information over the
graph topology. Wang & Leskovec (2020) used label propagation to improve GCNs.
Klicpera et al. (2019) incorporated personalized PageRank with GCNs. Liu et
al. (2021) introduced a $l_{1}$ norm-based graph smoothing term to enhance the
local smoothnesss adaptivity of GNNs. Hamilton et al. (2017); Zeng et al.
(2020) proposed sampling and aggregation frameworks to extent GCNs to
inductive learning settings. Another variant of GNNs is graph attention
networks (Velickovic et al., 2018; Thekumparampil et al., 2018; Abu-El-Haija
et al., 2018), which use attention mechanisms to adaptively learn aggregation
weights based on the nodes features. There are many other works on GNNs (Pei
et al., 2020) (Ying et al., 2018; Xinyi & Chen, 2019; Velickovic et al., 2019;
Zeng et al., 2020), we refer to Zhou et al. (2020); Battaglia et al. (2018);
Wu et al. (2021) for a comprehensive review. Most GNN models implicitly assume
that the labels of nodes and their neighbors should be the same or consistent,
while it does not hold for heterophilic graphs. Zhu et al. (2020) investigated
the issues of GNNs on heterophilic graphs and proposed to separately learn the
embeddings of ego-node and its neighborhood. Zhu et al. (2021) proposed a
framework to model the heterophily or homophily levels of graphs. Chien et al.
(2021) incorporated generalized PageRank with graph convolution to adapt GNNs
to heterophilic graphs.

There are also some works on the interpretability of GNNs proposed recently.
Li et al. (2018); Ying et al. (2019); Fu et al. (2020) showed that spectral
graph convolutions work as conducting Laplacian smoothing on the graph signals
and Wu et al. (2019); NT & Maehara (2019) demonstrated that GCN, SGC work as
low-pass filters. Gama et al. (2020) studied the stability properties of GNNs.
Xu et al. (2019); Oono & Suzuki (2020); Loukas (2020) studied the
expressiveness of GNNs. Verma & Zhang (2019); Garg et al. (2020) work on the
generalization and representation power of GNNs.

Graph based Semi-supervised Learning. Graph-based semi-supervised learning
works under the assumption that the labels of a node and its neighbors shall
be the same or consistent. Many methods have been proposed in the last decade,
such as Smola & Kondor (2003); Zhou et al. (2003); Belkin et al. (2004) use
Laplacian regularization techniques to force the labels of linked nodes to be
the same or consistent. Zhou & Schölkopf (2005) introduce discrete
regularization techniques to impose different regularizations on the node
features based on $p$-Laplacian. Lable propagation (Zhu et al., 2003)
recursively propagates the labels of labeled nodes over the graph topology and
use the convergence results to make predictions. To mention but a few, we
refer to Zhou & Schölkopf (2005); van Engelen & Hoos (2020) for a more
comprehensive review.

##  Appendix B Discussions and Future Work

In this section, we discuss the future work of pGNNs. Our theoretical results
and experimental results could lead to several potential extensions of pGNNs.

New Paradigm of Designing GNN Architectures. We bridge the gap between
discrete regularization framework, graph-based semi-supervised learning, and
GNNs, which provides a new paradigm of designing new GNN architectures.
Following the new paradigm, researchers could introduce more regularization
techniques, e.g., Laplacian regularization (Smola & Kondor, 2003; Belkin et
al., 2004), manifold regularization (Sindhwani et al., 2005; Belkin et al.,
2005; Niyogi, 2013), high-order regularization (Zhou & Belkin, 2011), Bayesian
regularization (Liu et al., 2014), entropy regularization (Grandvalet &
Bengio, 2004), and consider more explicit assumptions on graphs, e.g. the
homophily assumption, the low-density region assumption (i.e. the decision
boundary is likely to lie in a low data density region), manifold assumption
(i.e. the high dimensional data lies on a low-dimensional manifold), to
develop new graph convolutions or message passing schemes for graphs with
specific properties and generalize GNNs to a much broader range of graphs.
Moreover, the paradigm also enables us to explicitly study the behaviors of
the designed graph convolutions or message passing schemes from the theory of
regularization (Belkin & Niyogi, 2008; Niyogi, 2013; Slepcev & Thorpe, 2017).

Applications of pGNNs to learn on graphs with noisy topologies. The empirical
results (as shown in Fig. 3 and Tables 6 and 7) on graphs with noisy edges
show that pGNNs are very robust to noisy edges, which suggests the
applications of $p$-Laplacian message passing and pGNNs on the graph learning
scenarios where the graph topology could potentially be seriously intervened.

Integrating with existing GNN architectures. As shown in Table 9, the
experimental results on heterophilic benchmark datasets illustrate that
integrating GCN, JKNet with pGNNs can significantly improve their performance
on heterophilic graphs. It shows that pGNN could be used as a plug-and-play
component to be integrated into existing GNN architectures and improve their
performance on real-world applications.

Inductive learning for pGNNs. pGNNs are shown to be very effective for
inductive learning on PPI datasets as reported in Table 10. pGNNs even
outperforms GAT on PPI, while using much fewer parameters than GAT. It
suggests the promising extensions of pGNNs to inductive learning on graphs.

##  Appendix C Additional Theorems

###  C.1 Theorem 4 (Upper-Bounding Risk of pGNN)

######  Theorem 4 (Upper-Bounding Risks of pGNNs).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$ with
$N$ nodes, let ${\mathbf{X}}\in\mathbb{R}^{N\times c}$ be the node features
and ${\mathbf{y}}\in\mathbb{R}^{N}$ be the node labels and
${\mathbf{M}}^{(k)},\bm{\alpha}^{(k)},\bm{\beta}^{k},{\mathbf{F}}^{k}$ are
updated accordingly by Equations 12, 13 and 14 for $k=0,1,\dots,K-1$ and
${\mathbf{F}}^{(0)}={\mathbf{X}}$, $K\in{\mathbb{N}}$. Assume that
${\mathcal{G}}$ is $d$-regular and the ground-truth node features
${\mathbf{X}}^{*}={\mathbf{X}}+\bm{\epsilon}$, where
$\bm{\epsilon}\in\mathbb{R}^{N\times c}$ represents the noise in the node
features and there exists a $L$-Lipschitz function $\sigma:\mathbb{R}^{N\times
c}\rightarrow\mathbb{R}^{N}$ such that
$\sigma({\mathbf{X}}^{*})={\mathbf{y}}$. let
$\tilde{{\mathbf{y}}}^{(k+1)}=\bm{\alpha}^{(k)}\bm{{\mathbf{D}}}^{-1/2}{\mathbf{M}}^{(k)}\bm{{\mathbf{D}}}^{-1/2}\sigma({\mathbf{F}}^{(k)})+\bm{\beta}^{(k)}\sigma\left({\mathbf{F}}^{(0)}\right)$,
we have

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left|{y}_{i}-\tilde{{y}}_{i}\right|\leq{}$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|$
$\displaystyle\qquad+\frac{L}{N}\sum_{i=1}^{N}\alpha_{i,i}^{(K-1)}\left\|\Delta_{p}^{(K-1)}{\mathbf{F}}_{i,:}^{(K-1)}+\sum_{k=0}^{K-2}\prod_{l=k}^{K-2}\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(l)}}{d}\right)\Delta_{p}^{(k)}{\mathbf{X}}_{i,:}\right\|$
$\displaystyle\qquad+\frac{L}{N}\sum_{i=1}^{N}(1-\beta_{i,i}^{(K-1)})\left\|\bm{\epsilon}_{i,:}\right\|.$

Proof see Sec. D.4. Theorem 4 shows that the risk of pGNNs is upper-bounded by
the sum of three terms: The first term of the r.h.s in the above inequation
represents the risk of label prediction using only the original node features
${\mathbf{X}}$, the second term is the norm of $p$-Laplacian diffusion on the
node features ${\mathbf{X}}$, and the third term is the magnitude of the noise
in the node features. $\alpha_{i,i}$ and $\beta_{i,i}$ control the trade-off
between these three terms and they are related to the hyperparameter $\mu$ in
Eq. 10. The smaller $\mu$, the smaller $\beta_{i,i}$ and larger
$\alpha_{i,i}$, thus the more important of the $p$-Laplacian diffusion term
but also the more effect from the noise. Therefore, for graphs whose
topological information is not helpful for label prediction, we could impose
more weights on the first term by using a large $\mu$ so that pGNNs work more
like MLPs which simply learn on node features. While for graphs whose
topological information is helpful for label prediction, we could impose more
weights on the second term by using a small $\mu$ so that pGNNs can benefit
from $p$-Laplacian smoothing on node features.

In practice, to choose a proper value of $\mu$ one may first simply apply MLPs
on the node features to have a glance at the helpfulness of the node features.
If MLPs work very well, there is not much space for the graph’s topological
information to further improve the prediction performance and we may choose a
large $\mu$. Otherwise, there could be a large chance for the graph’s
topological information to further improve the performance and we should
choose a small $\mu$.

###  C.2 Theorem 5 ($p$-Orthogonal Theorem (Luo et al., 2010))

######  Theorem 5 ($p$-Orthogonal Theorem (Luo et al., 2010)).

If ${\mathbf{u}}^{(l)}$ and ${\mathbf{u}}^{(r)}$ are two eigenvectors of
$p$-Laplacian $\Delta_{p}$ associated with two different non-zero eigenvalues
$\lambda_{l}$ and $\lambda_{r}$, ${\mathbf{W}}$ is symmetric and $p\geq 1$,
then ${\mathbf{u}}^{(l)}$ and ${\mathbf{u}}^{(r)}$ are $p$-orthogonal up to
the second order Taylor expansion.

Theorem 5 implies that
$\phi_{p}({\mathbf{u}})^{(l)\top}\phi_{p}({\mathbf{u}}^{(r)})\approx 0$, for
all $l,r=0,\dots,N-1$ and $\lambda_{l}\neq\lambda_{r}$. Therefore, the space
spanned by the multiple eigenvectors of the graph $p$-Laplacian is
$p$-orthogonal.

###  C.3 Theorem 6 ($p$-Eigen-Decomposition of $\Delta_{p}$)

######  Theorem 6 ($p$-Eigen-Decomposition of $\Delta_{p}$).

Given the $p$-eigenvalues $\\{\lambda_{l}\in\mathbb{R}\\}_{l=0,1,\dots,N-1}$,
and the $p$-eigenvectors
$\\{{\mathbf{u}}^{(l)}\in\mathbb{R}^{N}\\}_{l=0,1,\dots,N-1}$ of $p$-Laplacian
$\Delta_{p}$ and
$\|{\mathbf{u}}^{(l)}\|_{p}=(\sum_{i=1}^{N}({u}_{i}^{(l)})^{p})^{1/p}=1$, let
${\mathbf{U}}$ be a matrix of $p$-eigenvectors with
${\mathbf{U}}=({\mathbf{u}}^{(0)},{\mathbf{u}}^{(1)},\dots,{\mathbf{u}}^{(N-1)})$
and ${\mathbf{\Lambda}}$ be a diagonal matrix with
${\mathbf{\Lambda}}=\mathrm{diag}(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1})$,
then the $p$-eigen-decomposition of $p$-Laplacian $\Delta_{p}$ is given by

$\displaystyle\Delta_{p}=\Phi_{p}({\mathbf{U}})\Lambda\Phi_{p}({\mathbf{U}})^{\top}.$

When $p=2$, it reduces to the standard eigen-decomposition of the Laplacian
matrix.

Proof see Sec. D.5.

###  C.4 Theorem 7 (Bounds of $p$-Eigenvalues)

######  Theorem 7 (Bounds of $p$-Eigenvalues).

Given a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}},{\mathbf{W}})$, if
${\mathcal{G}}$ is connected and $\lambda$ is a $p$-eigenvalue associated with
the $p$-eigenvector ${\mathbf{u}}$ of $\Delta_{p}$, let $N_{i}$ denotes the
number of edges connected to node $i$,
$N_{\text{min}}=\min\\{N_{i}\\}_{i=1,2,\dots,N}$, and
$k=\operatorname*{arg\,max}(\\{|{u}_{i}|/\sqrt{{D}_{i,i}}\\}_{i=1,2,\dots,N})$,
then

  1. 1.

for $p\geq 2$, $0\leq\lambda\leq 2^{p-1}$;

  2. 2.

for $1<p<2$, $0\leq\lambda\leq 2^{p-1}\sqrt{N_{k}}$;

  3. 3.

for $p=1$, $0\leq\lambda\leq\sqrt{N_{\text{min}}}$.

Proof see Sec. D.6.

##  Appendix D Proof of Theorems

###  D.1 Proof of Theorem 1

###### Proof.

Let $\bm{i}$ be the one-hot indicator vector whose $i$-th element is one and
the other elements are zero. Then, we can obtain the personalized PageRank on
node $i$, denoted as $\bm{\pi}_{\text{PPR}}(\bm{i})$, by using the recurrent
equation (Klicpera et al., 2019):

$\bm{\pi}_{\text{PPR}}^{(k+1)}(\bm{i})=\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\bm{\pi}_{\text{PPR}}^{(k)}(\bm{i})+\beta\bm{i},$

where $k$ is the iteration step, $0<\alpha<1$ and $\beta=(1-\alpha)$
represents the restart probability. Without loss of generality, suppose
$\bm{\pi}_{\text{PPR}}^{(0)}(\bm{i})=\bm{i}$. Then we have,

$\displaystyle\bm{\pi}_{\text{PPR}}^{(k)}(\bm{i})={}$
$\displaystyle\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\bm{\pi}_{\text{PPR}}^{(k-1)}(\bm{i})+\beta\bm{i}$
$\displaystyle={}$
$\displaystyle\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\left(\alpha\bm{D}^{-1/2}\bm{W}\bm{D}^{-1/2}\bm{\pi}_{\text{PPR}}^{(k-2)}(\bm{i})+\beta\bm{i}\right)+\beta\bm{i}$
$\displaystyle={}$
$\displaystyle\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{2}\bm{\pi}_{\text{PPR}}^{(t-2)}(\bm{i})+\beta\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\bm{i}+\beta\bm{i}$
$\displaystyle={}$
$\displaystyle\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{k}\bm{\pi}_{\text{PPR}}^{(0)}(\bm{i})+\beta\sum_{t=0}^{k-1}\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{t}\bm{i}$
$\displaystyle={}$
$\displaystyle\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{k}\bm{i}+\beta\sum_{t=0}^{k-1}\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{t}\bm{i}$

Since $0<\alpha<1$ and the eigenvalues of
${\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}$ in $[-1,1]$, we have

$\lim_{k\rightarrow\infty}\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{k}=0,$

and we also have

$\lim_{k\rightarrow\infty}\sum_{t=0}^{k-1}\left(\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{t}=\left({\mathbf{I}}_{N}-\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{-1}.$

Therefore,

$\displaystyle\bm{\pi}_{\text{PPR}}(\bm{i})=\lim_{k\rightarrow\infty}\bm{\pi}_{\text{PPR}}^{(k)}(\bm{i})={}$
$\displaystyle\beta\left({\mathbf{I}}_{N}-\alpha{\mathbf{D}}^{-1/2}{\mathbf{W}}{\mathbf{D}}^{-1/2}\right)^{-1}\bm{i}$
$\displaystyle={}$
$\displaystyle\beta\left(\alpha\Delta_{2}+(1-\alpha){\mathbf{I}}_{N}\right)^{-1}\bm{i}$
$\displaystyle={}$
$\displaystyle\mu(\Delta_{2}+\mu{\mathbf{I}}_{N})^{-1}\bm{i},$

where we let $\alpha=\frac{1}{1+\mu}$ and $\beta=\frac{\mu}{1+\mu}$, $\mu>0$.
Then the fully personalized PageRank matrix can be obtained by substituting
$\bm{i}$ with ${\mathbf{I}}_{N}$:

$\bm{\Pi}_{\text{PPR}}=\mu(\Delta_{2}+\mu{\mathbf{I}}_{N})^{-1}.$

∎

###  D.2 Proof of Theorem 2

###### Proof.

By the definition of ${\mathcal{L}}_{p}(f)$ in Eq. 8, we have for some
positive real value $\mu,\mu>0$

${\mathcal{L}}_{p}({\mathbf{F}})=\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}\right\|^{p}+\mu\sum_{i=1}^{N}\|{\mathbf{F}}_{i,:}-{\mathbf{X}}_{i,:}\|^{2}.$

and by Eq. 12,

${M}_{i,j}^{(k)}:={W}_{i,j}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|^{p-2}$

Then, we have

$\displaystyle\frac{\partial{\mathcal{L}}_{p}({\mathbf{F}}^{(k)})}{\partial{\mathbf{F}}_{i,:}^{(k)}}={}$
$\displaystyle
p\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|^{p-2}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right)+2\mu({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{X}}_{i,:})$
$\displaystyle={}$ $\displaystyle
p\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}{\mathbf{F}}_{i,:}^{(k)}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right)+2\mu({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{X}}_{i,:})$
$\displaystyle={}$ $\displaystyle
p\left(\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+\frac{2\mu}{p}\right){\mathbf{F}}_{i,:}^{(k)}-\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\frac{2\mu}{p}{\mathbf{X}}_{i,:}\right)\right)$
$\displaystyle={}$
$\displaystyle\frac{p}{\alpha_{i,i}^{(k)}}\left({\mathbf{F}}_{i,:}^{(k)}-\left(\alpha_{i,i}^{(k)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\beta_{i,i}^{(k)}{\mathbf{X}}_{i,:}\right)\right)$
$\displaystyle={}$
$\displaystyle\frac{p}{\alpha_{i,i}^{(k)}}\left({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{F}}_{i,:}^{(k+1)}\right),$

which indicates that

${\mathbf{F}}_{i,:}^{(k)}-{\mathbf{F}}_{i,:}^{(k+1)}=\frac{\alpha_{i,i}^{(k)}}{p}\cdot\frac{\partial{\mathcal{L}}_{p}({\mathbf{F}}^{(k)})}{\partial{\mathbf{F}}_{i,:}^{(k)}}.$

For all $i,j\in[N],{\mathbf{v}}\in\mathbb{R}^{1\times c}$, denote by

$\displaystyle\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}):={}$
$\displaystyle\frac{\partial{\mathcal{L}}_{p}({\mathbf{F}}^{(k)})}{\partial{\mathbf{F}}_{i,:}^{(k)}},$
$\displaystyle{M}_{i,j}^{\prime(k)}:={}$
$\displaystyle{W}_{i,j}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}({\mathbf{F}}_{i,:}^{(k)}+{\mathbf{v}})-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|^{p-2},$
$\displaystyle\alpha_{i,i}^{\prime(k)}:={}$ $\displaystyle
1\left/\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{{D}_{i,i}}+\frac{2\mu}{p}\right)\right.,$
$\displaystyle\beta_{i,i}^{\prime(k)}:={}$
$\displaystyle\frac{2\mu}{p}\alpha_{i,i}^{\prime(k)}$
$\displaystyle{\mathbf{F}}_{i,:}^{\prime(k+1)}:={}$
$\displaystyle\alpha_{i,i}^{\prime(k)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\beta_{i,i}^{\prime}{\mathbf{X}}_{i,:}.$

Then

$\displaystyle\left\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+{\mathbf{v}})-\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\right\|$
$\displaystyle={}$
$\displaystyle\left\|\frac{p}{\alpha_{i,i}^{\prime(k)}}\left({\mathbf{F}}_{i,:}^{(k)}+{\mathbf{v}}-{\mathbf{F}}_{i,:}^{\prime(k+1)}\right)-\frac{p}{\alpha_{i,i}^{(k)}}\left({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{F}}_{i,:}^{(k+1)}\right)\right\|$
$\displaystyle\leq{}$
$\displaystyle\frac{p}{\alpha_{i,i}^{\prime(k)}}\|{\mathbf{v}}\|+\left\|\frac{p}{\alpha_{i,i}^{\prime(k)}}\left({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{F}}_{i,:}^{\prime(k)}\right)-\frac{p}{\alpha_{i,i}^{(k)}}\left({\mathbf{F}}_{i,:}^{(k)}-{\mathbf{F}}_{i,:}^{(k+1)}\right)\right\|$
$\displaystyle={}$
$\displaystyle\frac{p}{\alpha_{i,i}^{\prime(k)}}\|{\mathbf{v}}\|+\left\|\left(\frac{p}{\alpha_{i,i}^{\prime(k)}}-\frac{p}{\alpha_{i,i}^{(k)}}\right){\mathbf{F}}_{i,:}^{(k)}-\frac{p}{\alpha_{i,i}^{\prime(k)}}{\mathbf{F}}_{i,:}^{\prime(k+1)}+\frac{p}{\alpha_{i,i}^{(k)}}{\mathbf{F}}_{i,:}^{(k+1)}\right\|$
$\displaystyle={}$
$\displaystyle\frac{p}{\alpha_{i,i}^{\prime(k)}}\|{\mathbf{v}}\|+p\left\|\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{{D}_{i,i}}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}\right){\mathbf{F}}_{i,:}^{(k)}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}-\frac{2\mu}{p}{\mathbf{X}}_{i,:}+\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\frac{2\mu}{p}{\mathbf{X}}_{i,:}\right\|$
$\displaystyle={}$
$\displaystyle\frac{p}{\alpha_{i,i}^{\prime(k)}}\|{\mathbf{v}}\|+p\left\|\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{{D}_{i,i}}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}\right){\mathbf{F}}_{i,:}^{(k)}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}+\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|$
$\displaystyle={}$
$\displaystyle\left(p\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+2\mu\right)\|{\mathbf{v}}\|+p\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}-{M}_{i,j}^{(k)}}{{D}_{i,i}}\|{\mathbf{v}}\|+p\left\|\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}-{M}_{i,j}^{(k)}}{{D}_{i,i}}{\mathbf{F}}_{i,:}^{(k)}-\sum_{j=1}^{N}\frac{{M}_{i,j}^{\prime(k)}-{M}_{i,j}^{(k)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|$
$\displaystyle={}$ $\displaystyle
p\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+\frac{2\mu}{p}+o\left(p,{\mathbf{v}},{\mathbf{X}},{\mathcal{G}}\right)\right)\|{\mathbf{v}}\|.$

Therefore, there exists some real positive value $\mu\in
o\left(p,{\mathbf{v}},{\mathbf{X}},{\mathcal{G}}\right)>0$ such that

$\left\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+{\mathbf{v}})-\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\right\|\leq
p\big{(}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+\frac{2\mu}{p}\big{)}\|{\mathbf{v}}\|=\frac{p}{\alpha_{i,i}^{(k)}}\|{\mathbf{v}}\|.$
(19)

Let $\bm{\gamma}=(\gamma_{1},\dots,\gamma_{N})^{\top}\in\mathbb{R}^{N}$ and
$\bm{\eta}\in\mathbb{R}^{N\times c}$. By Taylor’s theorem, we have:

$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+\gamma_{i}\bm{\eta}_{i,:})$
$\displaystyle={}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})+\gamma_{i}\int_{0}^{1}\langle\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+\epsilon\gamma_{i}\bm{\eta}_{i,:}),\bm{\eta}_{i,:}\rangle\mathrm{d}\epsilon$
$\displaystyle={}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})+\gamma_{i}\langle\bm{\eta}_{i,:},\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\rangle+\gamma_{i}\int_{0}^{1}\langle\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+\epsilon\gamma_{i}\bm{\eta}_{i,:})-\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}),\bm{\eta}_{i,:}\rangle\mathrm{d}\epsilon$
$\displaystyle\leq{}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})+\gamma_{i}\langle\bm{\eta}_{i,:},\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\rangle+\gamma_{i}\int_{0}^{1}\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}+\epsilon\gamma_{i}\bm{\eta}_{i,:})-\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\|\|\bm{\eta}_{i,:}\|\mathrm{d}\epsilon$
$\displaystyle\leq{}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})+\gamma_{i}\langle\bm{\eta}_{i,:},\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\rangle+\frac{p}{2\alpha_{i,i}^{(k)}}\gamma_{i}^{2}\|\bm{\eta}_{i,:}\|^{2}$

Let $\bm{\eta}=-\nabla{\mathcal{L}}_{p}({\mathbf{F}}^{(k)})$ and choose some
positive real value $\mu$ which depends on ${\mathbf{X}},{\mathcal{G}},p$ and
$p>1$, i.e. $\mu\in o\left(p,{\mathbf{X}},{\mathcal{G}}\right)$. By Eq. 19, we
have for all $i\in[N]$,

$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}-\gamma\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}))\leq{}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})-\langle\gamma_{i}\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)}),\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\rangle+\frac{p}{2\alpha_{i,i}^{(k)}}\gamma_{i}^{2}\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\|^{2}$
$\displaystyle={}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})-\frac{p}{2\alpha_{i,i}^{(k)}}\left(\frac{2\alpha_{i,i}^{(k)}\gamma_{i}}{p}-\gamma_{i}^{2}\right)\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:})\|^{2}$
$\displaystyle={}$
$\displaystyle{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})-\frac{p}{2\alpha_{i,i}^{(k)}}\left(\frac{\left(\alpha_{i,i}^{(k)}\right)^{2}}{p^{2}}-\left(\gamma_{i}-\frac{\alpha_{i,i}^{(k)}}{p}\right)^{2}\right)\|\partial{\mathcal{L}}_{p}({\mathbf{F}}_{i,:}^{(k)})\|^{2}.$

Then for all $i\in[N]$, when
$0\leq\gamma_{i}\leq\frac{2\alpha_{i,i}^{(k)}}{p}$, we have
${\mathcal{S}}_{p}({\mathbf{F}}_{i,:}^{(k)}-\gamma_{i}\partial{\mathcal{S}}_{p}({\mathbf{F}}_{i,:}^{(k)}))\leq{\mathcal{S}}_{p}({\mathbf{F}}_{i,:}^{(k)})$
and $\gamma_{i}=\frac{\alpha_{i,i}^{(k)}}{p}$ minimizes
${\mathcal{S}}_{p}({\mathbf{F}}_{i,:}^{(k)}-\gamma_{i}\partial{\mathcal{S}}_{p}({\mathbf{F}}_{i,:}^{(k)}))$.
Therefore,

${\mathcal{L}}_{p}({\mathbf{F}}^{(k+1)})={\mathcal{L}}_{p}({\mathbf{F}}^{(k)}-\frac{1}{p}\cdot\bm{\alpha}^{(k)}\nabla{\mathcal{L}}_{p}({\mathbf{F}}^{(k)}))\leq{\mathcal{L}}_{p}({\mathbf{F}}^{(k)}).$

∎

###  D.3 Proof of Theorem 3

###### Proof.

Without loss of generality, suppose ${\mathbf{F}}^{(0)}={\mathbf{X}}$. Denote
$\tilde{{\mathbf{M}}}^{(k)}={\mathbf{D}}^{-1/2}{\mathbf{M}}^{(k)}{\mathbf{D}}^{-1/2}$,
by Eq. 14, we have for $K\geq 2$,

$\displaystyle{\mathbf{F}}^{(K)}={}$
$\displaystyle\bm{\alpha}^{(K-1)}{\mathbf{D}}^{-1/2}{\mathbf{M}}^{(K-1)}{\mathbf{D}}^{-1/2}{\mathbf{F}}^{(K-1)}+\bm{\beta}^{(K-1)}{\mathbf{X}}$
$\displaystyle={}$
$\displaystyle\bm{\alpha}^{(K-1)}\tilde{{\mathbf{M}}}^{K-1}{\mathbf{F}}^{(K-1)}+\bm{\beta}^{(K-1)}{\mathbf{X}}$
$\displaystyle={}$
$\displaystyle\bm{\alpha}^{(K-1)}\tilde{{\mathbf{M}}}^{K-1}\left(\bm{\alpha}^{(K-2)}\tilde{{\mathbf{M}}}^{K-2}{\mathbf{F}}^{(K-2)}+\bm{\beta}^{(K-2)}{\mathbf{X}}\right)+\bm{\beta}^{(K-1)}{\mathbf{X}}$
$\displaystyle={}$
$\displaystyle\bm{\alpha}^{(K-1)}\bm{\alpha}^{(K-2)}\tilde{{\mathbf{M}}}^{K-1}\tilde{{\mathbf{M}}}^{K-2}{\mathbf{F}}^{(K-2)}+\bm{\alpha}^{(K-1)}\tilde{{\mathbf{M}}}^{K-1}\bm{\beta}^{(K-2)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}$
$\displaystyle={}$
$\displaystyle\left(\prod_{k=0}^{K-1}\bm{\alpha}^{(k)}\right)\left(\prod_{k=0}^{K-1}\tilde{{\mathbf{M}}}^{(k)}\right){\mathbf{F}}^{(0)}+\sum_{k=1}^{K-1}\left(\prod_{l=K-k}^{K-1}\bm{\alpha}^{(l)}\tilde{{\mathbf{M}}}^{(l)}\right)\bm{\beta}^{(K-1-k)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}$
$\displaystyle={}$
$\displaystyle\left(\prod_{k=0}^{K-1}\bm{\alpha}^{(k)}\right)\left(\prod_{k=0}^{K-1}\tilde{{\mathbf{M}}}^{(k)}\right){\mathbf{X}}+\sum_{k=1}^{K-1}\left(\prod_{l=K-k}^{K-1}\bm{\alpha}^{(l)}\tilde{{\mathbf{M}}}^{(l)}\right)\bm{\beta}^{(K-1-k)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}.$
(20)

Recall Equations 12 and 13, we have

$\tilde{{M}}_{i,j}^{(k)}=\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}^{(k)}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}^{(k)}\right\|^{p-2},\text{
for all }i,j=1,2,\dots,N,$

and

$\alpha_{i,i}^{(k)}=\frac{1}{{\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}+\frac{2\mu}{p}}},\text{
for all }i=1,2,\dots,N.$

Note that the eigenvalues of $\tilde{{\mathbf{M}}}$ are not infinity and
$0<\alpha_{i,i}<1$ for all $i=1,\dots,N$. Then we have

$\lim_{K\rightarrow\infty}\prod_{k=0}^{K-1}\bm{\alpha}^{(k)}=0,$

and

$\lim_{K\rightarrow\infty}\left(\prod_{k=0}^{K-1}\bm{\alpha}^{(k)}\right)\left(\prod_{k=0}^{K-1}\tilde{{\mathbf{M}}}^{(k)}\right)=0.$

Therefore,

$\lim_{K\rightarrow\infty}{\mathbf{F}}^{(K)}=\lim_{K\rightarrow\infty}\left(\sum_{k=1}^{K-1}\left(\prod_{l=K-k}^{K-1}\bm{\alpha}^{(l)}\tilde{{\mathbf{M}}}^{(l)}\right)\bm{\beta}^{(K-1-k)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}\right).$
(21)

By Equations 6 and 12, we have

$\displaystyle\Delta_{p}f(i)={}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\|(\nabla
f)([j,i])\|^{p-2}f(i)-\sum_{j=1}^{N}\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\|(\nabla
f)([j,i])\|^{p-2}f(j)$ $\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}f(i)-\sum_{j=1}^{N}\frac{{M}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}f(j).$
(22)

By Eq. 13, we have

$\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}=\frac{1}{\alpha_{i,i}^{(k)}}-\frac{2\mu}{p}.$
(23)

Equations 22 and 23 show that

$\Delta_{p}^{(k)}=\left(\left(\bm{\alpha}^{(k)}\right)^{-1}-\frac{2\mu}{p}{\mathbf{I}}_{N}\right)-\tilde{{\mathbf{M}}}^{(k)},$
(24)

which indicates

$\bm{\alpha}^{(k)}\tilde{{\mathbf{M}}}^{(k)}={\mathbf{I}}_{N}-\frac{2\mu}{p}\bm{\alpha}^{(k)}-\bm{\alpha}^{(k)}\Delta_{p}^{(k)}.$
(25)

Eq. 25 shows that $\bm{\alpha}^{(k)}\tilde{{\mathbf{M}}}^{(k)}$ is linear
w.r.t $\Delta_{p}$ and therefore can be expressed by a linear combination in
terms of $\Delta_{p}$:

$\bm{\alpha}^{(k)}\tilde{{\mathbf{M}}}^{(k)}=\bm{\theta}^{\prime(k)}\Delta_{p},$
(26)

where
$\bm{\theta}^{\prime}=\mathrm{diag}(\theta_{0}^{\prime},\theta_{1}^{\prime},\dots,\theta_{N-1}^{\prime})$
are the parameters. Therefore, we have

$\displaystyle\lim_{K\rightarrow\infty}{\mathbf{F}}^{(K)}={}$
$\displaystyle\lim_{K\rightarrow\infty}\left(\sum_{k=1}^{K-1}\left(\prod_{l=K-k}^{K-1}\bm{\alpha}^{(l)}\tilde{{\mathbf{M}}}^{(l)}\right)\bm{\beta}^{(K-1-k)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}\right)$
$\displaystyle={}$
$\displaystyle\lim_{K\rightarrow\infty}\left(\sum_{k=1}^{K-1}\left(\prod_{l=K-k}^{K-1}\bm{\theta}^{\prime(l)}\Delta_{p}\right)\bm{\beta}^{(K-1-k)}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}\right)$
$\displaystyle={}$
$\displaystyle\lim_{K\rightarrow\infty}\left(\sum_{k=1}^{K-1}\bm{\beta}^{(K-1-k)}\left(\prod_{l=K-k}^{K-1}\bm{\theta}^{\prime(l)}\right)\Delta_{p}^{k}{\mathbf{X}}+\bm{\beta}^{(K-1)}{\mathbf{X}}\right)$
$\displaystyle={}$
$\displaystyle\lim_{K\rightarrow\infty}\sum_{k=0}^{K-1}\bm{\theta}^{\prime\prime(k)}\Delta_{p}^{k}{\mathbf{X}},$

where
$\bm{\theta}^{\prime\prime(k)}=\mathrm{diag}(\theta_{1}^{\prime\prime(k)},\theta_{2}^{\prime\prime(k)},\dots,\theta_{N}^{\prime\prime(k)})$
defined as $\bm{\theta}^{\prime\prime(0)}=\bm{\beta}^{(K-1)}$ and

$\theta_{i}^{\prime\prime(k)}=\beta_{i,i}^{(K-1-k)}\prod_{l=K-k}^{K-1}\theta_{i}^{\prime(l)},\text{
for }k=1,2,\dots,K-1.$

Let $\bm{\theta}=(\theta_{0},\theta_{1},\dots,\theta_{K-1})$ defined as
$\theta_{k}=\sum_{i=1}^{N}\theta_{i}^{\prime\prime(k)}$ for all
$k=0,1,\dots,K-1$, then

$\displaystyle\lim_{K\rightarrow\infty}{\mathbf{F}}^{(K)}={}$
$\displaystyle\lim_{K\rightarrow\infty}\left(\sum_{k=1}^{K-1}\theta_{k}\Delta_{p}^{k}{\mathbf{X}}+\theta_{0}{\mathbf{X}}\right)$
$\displaystyle={}$
$\displaystyle\lim_{K\rightarrow\infty}\sum_{k=0}^{K-1}\theta_{k}\Delta_{p}^{k}{\mathbf{X}}.$

Therefore complete the proof. ∎

###  D.4 Proof of Theorem 4

###### Proof.

The first-order Taylor expansion with Peano’s form of remainder for $\sigma$
at ${\mathbf{X}}_{i,:}^{*}$ is given by:

$\sigma({\mathbf{F}}_{j,:}^{(K-1)})=\sigma({\mathbf{X}}_{i,:}^{*})+\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}^{*}\right)^{\top}+o(\|{\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}^{*}\|).$

Note that in general the output non-linear layer $\sigma(\cdot)$ is simple.
Here we assume that it can be well approximated by the first-order Taylor
expansion and we can ignore the Peano’s form of remainder. For all
$i=1,\dots,N$, ${D}_{i,i}={D}_{j,j}=d$, we have
$\alpha_{i,i}\sum_{j=1}^{N}{D}_{i,i}^{-1/2}{M}_{i,j}^{(K-1)}{D}_{j,j}^{-1/2}+\beta_{i,i}^{(K-1)}=1$.
Then

$\displaystyle\left|{y}_{i}-\tilde{{y}}_{i}^{(K)}\right|$ $\displaystyle={}$
$\displaystyle\left|{y}_{i}-\alpha_{i,i}^{(K-1)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\sigma\left({\mathbf{F}}_{j,:}^{K-1}\right)-\beta_{i,i}^{(K-1)}\sigma\left({\mathbf{X}}_{i,:}\right)\right|$
$\displaystyle={}$
$\displaystyle\left|{y}_{i}-\beta_{i,i}^{(K-1)}\sigma\left({\mathbf{X}}_{i,:}\right)-\alpha_{i,i}^{(K-1)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\left(\sigma\left({\mathbf{X}}_{i,:}^{*}\right)+\frac{\partial\sigma\left({\mathbf{X}}_{i,:}^{*}\right)}{\partial{\mathbf{X}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}^{*}\right)^{\top}\right)\right|$
$\displaystyle={}$
$\displaystyle\left|{y}_{i}-\alpha_{i,i}^{(K-1)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}{y}_{i}-\beta_{i,i}^{(K-1)}\sigma({\mathbf{X}}_{i,:})-\alpha_{i,i}^{(K-1)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\left(\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}^{*}\right)^{\top}\right)\right|$
$\displaystyle={}$
$\displaystyle\left|\beta_{i,i}^{(K-1)}({y}_{i}-\sigma({\mathbf{X}}_{i,:}))-\alpha_{i,i}^{(K-1)}\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\left(\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}-\bm{\epsilon}_{i,:}\right)^{\top}\right)\right|$
$\displaystyle\leq{}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}\left|\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}\right)^{\top}\right|+(1-\beta_{i,i}^{(K-1)})\left|\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\bm{\epsilon}_{i,:}^{\top}\right|$
$\displaystyle\leq{}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}\left\|\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\right\|\left\|\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}\right)\right\|+(1-\beta_{i,i}^{(K-1)})\left\|\frac{\partial\sigma({\mathbf{X}}_{i,:}^{*})}{\partial{\mathbf{X}}}\right\|\left\|\bm{\epsilon}_{i,:}\right\|$
$\displaystyle\leq{}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}L\left\|\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{d}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{X}}_{i,:}\right)\right\|+(1-\beta_{i,i}^{(K-1)})L\left\|\bm{\epsilon}_{i,:}\right\|$
$\displaystyle={}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}L\left\|\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{d}\left({\mathbf{F}}_{j,:}^{(K-1)}-{\mathbf{F}}_{i,:}^{(K-1)}+{\mathbf{F}}_{i,:}^{(K-1)}-{\mathbf{X}}_{i,:}\right)\right\|$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+(1-\beta_{i,i}^{(K-1)})L\left\|\bm{\epsilon}_{i,:}\right\|$
$\displaystyle={}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}L\left\|\Delta_{p}{\mathbf{F}}_{i,:}^{(K-1)}+\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-1)}}{d}\cdot\sum_{j=1}^{N}\frac{{M}_{i,j}^{(K-2)}}{d}\left({\mathbf{F}}_{j,:}^{(K-2)}-{\mathbf{X}}_{i,:}\right)\right\|$
$\displaystyle\qquad\qquad\qquad\qquad\qquad+(1-\beta_{i,i}^{(K-1)})L\left\|\bm{\epsilon}_{i,:}\right\|$
$\displaystyle={}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}L\left\|\Delta_{p}^{(K-1)}{\mathbf{F}}_{i,:}^{(K-1)}+\sum_{k=0}^{K-2}\prod_{l=k}^{K-2}\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(l)}}{d}\right)\Delta_{p}^{(k)}{\mathbf{X}}_{i,:}\right\|+(1-\beta_{i,i}^{(K-1)})L\left\|\bm{\epsilon}_{i,:}\right\|$
$\displaystyle={}$
$\displaystyle\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|+\alpha_{i,i}^{(K-1)}L\left\|\Delta_{p}^{(K-1)}{\mathbf{F}}_{i,:}^{(K-1)}+\sum_{k=0}^{K-2}\prod_{l=k}^{K-2}\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(l)}}{d}\right)\Delta_{p}^{(k)}{\mathbf{X}}_{i,:}\right\|+(1-\beta_{i,i}^{(K-1)})L\left\|\bm{\epsilon}_{i,:}\right\|$

Therefore,

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\left|{y}_{i}-\tilde{{y}}_{i}\right|\leq{}$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\beta_{i,i}^{(K-1)}\left|{y}_{i}-\sigma({\mathbf{X}}_{i,:})\right|$
$\displaystyle\qquad+\frac{L}{N}\sum_{i=1}^{N}\alpha_{i,i}^{(K-1)}\left\|\Delta_{p}^{(K-1)}{\mathbf{F}}_{i,:}^{(K-1)}+\sum_{k=0}^{K-2}\prod_{l=k}^{K-2}\left(\sum_{j=1}^{N}\frac{{M}_{i,j}^{(l)}}{d}\right)\Delta_{p}^{(k)}{\mathbf{X}}_{i,:}\right\|$
$\displaystyle\qquad+\frac{L}{N}\sum_{i=1}^{N}(1-\beta_{i,i}^{(K-1)})\left\|\bm{\epsilon}_{i,:}\right\|$

∎

###  D.5 Proof of Theorem 6

###### Proof.

Note that

$\phi_{p}({\mathbf{u}})^{\top}{\mathbf{u}}=\sum_{i=1}^{N}\phi_{p}({u}_{i}){u}_{i}=\sum_{i=1}^{N}\|{u}_{i}\|^{p-2}{u}_{i}^{2}=\sum_{i=1}^{N}\|{u}_{i}\|^{p}=\sum_{i=1}^{N}|{u}_{i}|^{p}=\|{\mathbf{u}}\|_{p}^{p}=1,$

then we have

$\Delta_{p}{\mathbf{U}}=\Phi_{p}({\mathbf{U}}){\mathbf{\Lambda}}=\Phi_{p}({\mathbf{U}}){\mathbf{\Lambda}}\Phi({\mathbf{U}})^{\top}{\mathbf{U}}.$

Therefore,
$\Delta_{p}=\Phi_{p}({\mathbf{U}}){\mathbf{\Lambda}}\Phi_{p}({\mathbf{U}})^{\top}$.

When $p=2$, by $\Phi_{2}({\mathbf{U}})={\mathbf{U}}$, we get
$\Delta_{2}=\Phi_{2}({\mathbf{U}}){\mathbf{\Lambda}}\Phi_{2}({\mathbf{U}})^{\top}={\mathbf{U}}{\mathbf{\Lambda}}{\mathbf{U}}^{\top}$.
∎

###  D.6 Proof of Theorem 7

###### Proof.

By the definition of graph $p$-Laplacian, we have for all $i=1,2,\dots,N$,

$\left(\Delta_{p}{\mathbf{u}}\right)_{i}=\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right\|^{p-2}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right)=\lambda\phi_{p}({u}_{i}).$

Then, for all $i=1,2,\dots,N$,

$\displaystyle\lambda={}$
$\displaystyle\frac{1}{\phi_{p}({u}_{i})}\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right\|^{p-2}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right)$
$\displaystyle={}$
$\displaystyle\frac{1}{\|{u}_{i}\|^{p-2}{u}_{i}}\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right\|^{p-2}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right)$
$\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\frac{\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right\|^{p-2}}{\|{u}_{i}\|^{p-2}}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}\frac{{u}_{j}}{{u}_{i}}\right)$
$\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{i,j}}{{D}_{i,i}}-\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\frac{{u}_{j}}{{u}_{i}}\right)\left(\frac{\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}\right\|}{\|{u}_{i}\|}\right)^{p-2}$
$\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{i,j}}{{D}_{i,i}}-\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\frac{{u}_{j}}{{u}_{i}}\right)\left\|\frac{\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{u}_{i}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{u}_{j}}{{u}_{i}}\right\|^{p-2}$
$\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{i,j}}{{D}_{i,i}}-\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\frac{{u}_{j}}{{u}_{i}}\right)\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}\frac{{u}_{j}}{{u}_{i}}\right\|^{p-2}$

Let $l=\operatorname*{arg\,max}\\{\|u_{i}\|\\}_{i=1,2,\dots,N}$, the above
equation holds for all $i=1,2,\dots,N$, then

$\displaystyle\lambda={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{l,j}}{{D}_{l,l}}-\frac{{W}_{l,j}}{\sqrt{{D}_{l,l}{D}_{j,j}}}\frac{{u}_{j}}{{u}_{l}}\right)\left\|\sqrt{\frac{{W}_{l,j}}{{D}_{l,l}}}-\sqrt{\frac{{W}_{l,j}}{{D}_{j,j}}}\frac{{u}_{j}}{{u}_{l}}\right\|^{p-2}$
$\displaystyle\geq{}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{l,j}}{{D}_{l,l}}-\frac{{W}_{l,j}}{\sqrt{{D}_{l,l}{D}_{j,j}}}\left|\frac{{u}_{j}}{{u}_{l}}\right|\right)\left\|\sqrt{\frac{{W}_{l,j}}{{D}_{l,l}}}-\sqrt{\frac{{W}_{l,j}}{{D}_{j,j}}}\frac{{u}_{j}}{{u}_{l}}\right\|^{p-2}$
$\displaystyle\geq{}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{l,j}}{{D}_{l,l}}-\frac{{W}_{l,j}}{\sqrt{{D}_{l,l}{D}_{j,j}}}\right)\left\|\sqrt{\frac{{W}_{l,j}}{{D}_{l,l}}}-\sqrt{\frac{{W}_{l,j}}{{D}_{j,j}}}\frac{{u}_{j}}{{u}_{l}}\right\|^{p-2}$
$\displaystyle\geq{}$ $\displaystyle 0.$

When $p=1$,

$\displaystyle\lambda={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{W_{i,j}}{D_{i,i}}-\frac{W_{i,j}}{\sqrt{D_{i,i}D_{j,j}}}\frac{u_{j}}{u_{i}}\right)\left\|\sqrt{\frac{W_{i,j}}{D_{i,i}}}-\sqrt{\frac{W_{i,j}}{D_{j,j}}}\frac{u_{j}}{u_{i}}\right\|^{-1}$
$\displaystyle\leq{}$
$\displaystyle\sum_{j=1}^{N}\sqrt{\frac{W_{i,j}}{D_{i,i}}}\leq\sqrt{N_{i}\sum_{j=1}^{N}\frac{W_{i,j}}{D_{i,i}}}=\sqrt{N_{i}},$

where the last inequality holds by using the Cauchy-Schwarz inequality. The
above inequality holds for all $i=1,2,\dots,N$, therefore,

$\lambda\leq\sum_{j=1}^{N}\sqrt{\frac{W_{i,j}}{D_{i,i}}}\leq\sqrt{N_{min}}.$

When $p>1$, we have for $i=1,2,\dots,N$,

$\displaystyle\lambda={}$
$\displaystyle\sum_{j=1}^{N}\left(\frac{{W}_{i,j}}{{D}_{i,i}}-\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\frac{{u}_{j}}{u_{i}}\right)\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}\frac{{u}_{j}}{u_{i}}\right\|^{p-2}$
$\displaystyle\leq{}$
$\displaystyle\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}\frac{{u}_{j}}{u_{i}}\right\|^{p-1}$
$\displaystyle\leq{}$
$\displaystyle\sum_{j=1}^{N}\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}+\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}\left|\frac{{u}_{j}}{u_{i}}\right|\right\|^{p-1}$
$\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}\right)^{p}\left\|1+\sqrt{\frac{{D}_{i,i}}{{D}_{j,j}}}\left|\frac{{u}_{j}}{u_{i}}\right|\right\|^{p-1}.$

Without loss of generality, let
$k=\operatorname*{arg\,max}(\\{|{u}_{i}|/\sqrt{{D}_{i,i}}\\}_{i=1,2,\dots,N})$.
Because the above inequality holds for all $i=1,2,\dots,N$, then we have

$\displaystyle\lambda\leq{}$
$\displaystyle\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\right)^{p}\left(1+\sqrt{\frac{{D}_{k,k}}{{D}_{j,j}}}\left|\frac{{u}_{j}}{u_{k}}\right|\right)^{p-1}$
$\displaystyle\leq{}$ $\displaystyle
2^{p-1}\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\right)^{p}.$

For $p\geq 2$,

$\lambda\leq
2^{p-1}\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\right)^{p}\leq
2^{p-1}\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\right)^{2}=2^{p-1}.$

For $1<p<2$,

$\lambda\leq
2^{p-1}\sum_{j=1}^{N}\left(\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\right)^{p}\leq
2^{p-1}\sum_{j=1}^{N}\sqrt{\frac{{W}_{k,j}}{{D}_{k,k}}}\leq
2^{p-1}\sqrt{N_{k}\sum_{j=1}^{N}\frac{W_{k,j}}{D_{k,k}}}=2^{p-1}\sqrt{N_{k}}.$

∎

###  D.7 Proof of Proposition 1

###### Proof.

We proof Proposition 1 based on the bounds of $p$-eigenvalues as demonstrated
in Theorem 7.

By Eq. 6 and Eq. 12, we have

$\displaystyle\Delta_{p}f(i)={}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\|(\nabla
f)([j,i])\|^{p-2}f(i)-\sum_{j=1}^{N}\frac{{W}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}\|(\nabla
f)([j,i])\|^{p-2}f(j)$ $\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}f(i)-\sum_{j=1}^{N}\frac{{M}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}f(j).$
(27)

By Eq. 13, we have

$\sum_{j=1}^{N}\frac{{M}_{i,j}^{(k)}}{{D}_{i,i}}=\frac{1}{\alpha_{i,i}^{(k)}}-\frac{2\mu}{p}.$
(28)

Equations 27 and 28 show that

$\Delta_{p}^{(k)}=\left(\left(\bm{\alpha}^{(k)}\right)^{-1}-\frac{2\mu}{p}{\mathbf{I}}_{N}\right)-{\mathbf{D}}^{-1/2}{\mathbf{M}}^{(k)}{\mathbf{D}}^{-1/2},$
(29)

which indicates

$\bm{\alpha}^{(k)}{\mathbf{D}}^{-1/2}{\mathbf{M}}^{(k)}{\mathbf{D}}^{-1/2}={\mathbf{I}}_{N}-\frac{2\mu}{p}\bm{\alpha}^{(k)}-\bm{\alpha}^{(k)}\Delta_{p}^{(k)}.$
(30)

For $i=1,2,\dots,N$, let
$\tilde{\bm{\alpha}}:=(\tilde{\alpha}_{1},\dots,\tilde{\alpha}_{N})$,
$\tilde{\alpha}_{i}:=1/\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}$, then

$\displaystyle\alpha_{i,i}\sum_{j=1}^{N}\frac{{M}_{i,j}}{\sqrt{{D}_{i,i}{D}_{j,j}}}={}$
$\displaystyle(1-\frac{2\mu}{p}\alpha_{i,i})-\alpha_{i,i}\lambda_{i}$
$\displaystyle={}$
$\displaystyle\frac{\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}}{\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}+\frac{2\mu}{p}}\left(1-\frac{1}{\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}}\lambda_{i}\right)$
$\displaystyle={}$
$\displaystyle\frac{1}{1+\frac{2\mu\tilde{\alpha}_{i}}{p}}\left(1-\tilde{\alpha}_{i}\lambda_{i}\right),$
(31)

Recall the Eq. 12 that

${M}_{i,j}={W}_{i,j}\left\|\sqrt{\frac{{W}_{i,j}}{{D}_{i,i}}}{\mathbf{F}}_{i,:}-\sqrt{\frac{{W}_{i,j}}{{D}_{j,j}}}{\mathbf{F}}_{j,:}\right\|^{p-2}={W}_{i,j}\|(\nabla
f)([i,j])\|^{p-2},$

  1. 1.

When $p=2$, for all $i=1,\dots,N$, $\tilde{\alpha}_{i}=1$ and
$0\leq\lambda_{i-1}\leq 2$, $g_{2}(\lambda_{i-1})$ works as low-high-pass
filters.

  2. 2.

When $p>2$, by Theorem 7 we have for all $i=1,\dots,N$,
$0\leq\lambda_{i-1}\leq 2^{p-1}$. If $0\leq\tilde{\alpha}_{i}\leq 2^{1-p}$,
then $0\leq 1-\tilde{\alpha}_{i}\lambda_{i}\leq 1$, which indicates that
$g_{p}(\lambda_{i-1})$ works as a low-pass filter; If
$\tilde{\alpha}_{i}>2^{1-p}$, then $g_{p}(\lambda_{i-1})$ works as low-high-
pass filters. Since

$\displaystyle\sum_{j=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}={}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}\|(\nabla
f)([i,j])\|^{p-2}}{{D}_{i,i}}$ $\displaystyle\leq{}$
$\displaystyle\sqrt{\sum_{j=1}^{N}\left(\frac{{W}_{i,j}}{{D}_{i,i}}\right)^{2}\sum_{j=1}^{N}\|(\nabla
f)([i,j])\|^{2(p-2)}}$ $\displaystyle\leq{}$
$\displaystyle\sqrt{\sum_{j=1}^{N}\|(\nabla f)([i,j])\|^{2(p-2)}}$
$\displaystyle\leq{}$ $\displaystyle\|\nabla f(i)\|^{p-2},$

which indicates that $\tilde{\alpha}_{i}\geq\|\nabla f(i)\|^{2-p}$.
$0\leq\tilde{\alpha}_{i}\leq 2^{1-p}$ directly implies that $0\leq\|\nabla
f(i)\|^{2-p}\leq 2^{1-p}$, i.e. $\|\nabla f(i)\|\geq 2^{(p-1)/(p-2)}$ and when
$\|\nabla f(i)\|^{2-p}\geq 2^{1-p}$, i.e. $\|\nabla f(i)\|\leq
2^{(p-1)/(p-2)}$, $\tilde{\alpha}_{i}\geq 2^{1-p}$ always holds. Therefore, if
$\|\nabla f(i)\|\leq 2^{(p-1)/(p-2)}$, $g_{p}(\lambda_{i-1})$ works as low-
high-pass filters on node $i$; If $g_{p}(\lambda_{i-1})$ works as a low-pass
filter, $\|\nabla f(i)\|\geq 2^{(p-1)/(p-2)}$.

  3. 3.

When $1\leq p<2$, by Theorem 7 we have for all $i=1,\dots,N$,
$0\leq\lambda_{i-1}\leq 2^{p-1}\sqrt{N_{k}}$. If $0\leq\tilde{\alpha}_{i}\leq
2^{1-p}/\sqrt{N_{k}}$, $0\leq 1-\tilde{\alpha}_{i}\lambda_{i}\leq 1$, which
indicates that $g_{p}(\lambda_{i-1})$ work as low-pass filters; If
$\tilde{\alpha}_{i}\geq 2^{1-p}/\sqrt{N_{k}}$, $g_{p}(\lambda_{i-1})$ work as
low-high-pass filters. By

$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\frac{1}{\left\|(\nabla
f)([i,j])\right\|^{p-2}}\left/\frac{1}{\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left\|(\nabla
f)([i,j])\right\|^{p-2}}\right.$ $\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\frac{1}{\left\|(\nabla
f)([i,j])\right\|^{p-2}}\cdot\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left\|(\nabla
f)([i,j])\right\|^{p-2}$ $\displaystyle\geq{}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left(\frac{1}{\left\|(\nabla
f)([i,j])\right\|^{p-2}}\cdot\left\|(\nabla f)([i,j])\right\|^{p-2}\right)$
$\displaystyle={}$ $\displaystyle 1,$

we have

$\displaystyle\tilde{\alpha}_{i}={}$
$\displaystyle\frac{1}{\sum_{i=1}^{N}\frac{{M}_{i,j}}{{D}_{i,i}}}=\frac{1}{\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left\|(\nabla
f)([i,j])\right\|^{p-2}}$ $\displaystyle\leq{}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\frac{1}{\left\|(\nabla
f)([i,j])\right\|^{p-2}}$ $\displaystyle={}$
$\displaystyle\sum_{j=1}^{N}\frac{{W}_{i,j}}{{D}_{i,i}}\left\|(\nabla
f)([i,j])\right\|^{2-p}$ $\displaystyle\leq{}$ $\displaystyle\|\nabla
f(i)\|^{2-p}.$

$\tilde{\alpha}_{i}\geq 2^{1-p}/\sqrt{N_{k}}$ directly implies that $\|\nabla
f(i)\|^{2-p}\geq 2^{1-p}/\sqrt{N_{k}}$, i.e. $\|\nabla f(i)\|\geq
2(2\sqrt{N_{k}})^{1/(p-2)}$ and when $0\leq\|\nabla f(i)\|^{2-p}\leq
2^{1-p}/\sqrt{N_{k}}$, i.e. $0\leq\|\nabla f(i)\|\leq
2(2\sqrt{N_{k}})^{1/(p-2)}$, $0\leq\tilde{\alpha}_{i}\leq
2^{1-p}/\sqrt{N_{k}}$ always holds. Therefore, if $0\leq\|\nabla f(i)\|\leq
2(2\sqrt{N_{k}})^{1/(p-2)}$, $g_{p}(\lambda_{i-1})$ work as low-pass filters;
If $g_{p}(\lambda_{i-1})$ work as low-high-pass filters, $\|\nabla f(i)\|\geq
2\left(2\sqrt{N_{k}}\right)^{1/(p-2)}$.

Specifically, when $p=1$, by Theorem 7 we have for all $i=1,\dots,N$,
$0\leq\lambda_{i-1}\leq 2^{p-1}\sqrt{N_{\text{min}}}$. Following the same
derivation above we attain if $0\leq\|\nabla f(i)\|\leq
2(2\sqrt{N_{\text{min}}})^{1/(p-2)}$, $g_{p}(\lambda_{i-1})$ work as low-pass
filters; If $g_{p}(\lambda_{i-1})$ work as both low-high-pass filters,
$\|\nabla f(i)\|\geq 2\left(2\sqrt{N_{\text{min}}}\right)^{1/(p-2)}$.

∎

##  Appendix E Dataset Statistics and Hyperparameters

###  E.1 Illustration of Graph Gradient and Graph Divergence

(a) Graph gradient.

(b) Graph divergence.

Figure 4: A tiny example of illustration of graph gradient and graph
divergence. Best view in colors.

###  E.2 Dataset Statistics

Table 2 summarizes the dataset statistics and the levels of homophily
${\mathcal{H}}({\mathcal{G}})$ of all benchmark datasets. Note that the
homophily scores here is different with the scores reported by Chien et al.
(2021). There is a bug in their code when computing the homophily scores
(doing division with torch integers) which caused their homophily scores to be
smaller. Additionally, We directly used the data from Pytorch Geometric
library (Fey & Lenssen, 2019) where they did not transform Chameleon and
Squirrel to undirected graphs, which is different from Chien et al. (2021)
where they did so.

Table 2: Statistics of datasets. Dataset | #Class | #Feature | #Node | #Edge | Training | Validation | Testing | ${\mathcal{H}}({\mathcal{G}})$
---|---|---|---|---|---|---|---|---
Cora | 7 | 1433 | 2708 | 5278 | 2.5% | 2.5% | 95% | 0.825
CiteSeer | 6 | 3703 | 3327 | 4552 | 2.5% | 2.5% | 95% | 0.717
PubMed | 3 | 500 | 19717 | 44324 | 2.5% | 2.5% | 95% | 0.792
Computers | 10 | 767 | 13381 | 245778 | 2.5% | 2.5% | 95% | 0.802
Photo | 8 | 745 | 7487 | 119043 | 2.5% | 2.5% | 95% | 0.849
CS | 15 | 6805 | 18333 | 81894 | 2.5% | 2.5% | 95% | 0.832
Physics | 5 | 8415 | 34493 | 247962 | 2.5% | 2.5% | 95% | 0.915
Chameleon | 5 | 2325 | 2277 | 31371 | 60% | 20% | 20% | 0.247
Squirrel | 5 | 2089 | 5201 | 198353 | 60% | 20% | 20% | 0.216
Actor | 5 | 932 | 7600 | 26659 | 60% | 20% | 20% | 0.221
Wisconsin | 5 | 251 | 499 | 1703 | 60% | 20% | 20% | 0.150
Texas | 5 | 1703 | 183 | 279 | 60% | 20% | 20% | 0.097
Cornell | 5 | 1703 | 183 | 277 | 60% | 20% | 20% | 0.386

###  E.3 Hyperparameter Settings

We set the number of layers as 2, the maximum number of epochs as 1000, the
number for early stopping as 200, the weight decay as 0 or 0.0005 for all
models. The other hyperparameters for each model are listed as below:

  * •

1.0GNN, 1.5GNN, 2.0GNN, 2.5GNN:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01,0.05\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

    * –

$\mu$: $\\{0.01,0.1,0.2,1,10\\}$

    * –

$K$: $4,6,8$

  * •

MLP:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

  * •

GCN:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

  * •

SGC:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.2,0.01\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

    * –

$K$: $2$

  * •

GAT:

    * –

Number of hidden units: $8$

    * –

Number of attention heads: $8$

    * –

Learning rate: $\\{0.001,0.005\\}$

    * –

Dropout rate: $\\{0,0.6\\}$

  * •

JKNet:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

    * –

$K$: $10$

    * –

$\alpha$: $\\{0.1,0.5,0.7,1\\}$

    * –

The number of GCN based layers: $2$

    * –

The layer aggregation: LSTM with $16$ channels and $4$ layers

  * •

APPNP:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

    * –

$K$: $10$

    * –

$\alpha$: $\\{0.1,0.5,0.7,1\\}$

  * •

GPRGNN:

    * –

Number of hidden units: $16$

    * –

Learning rate: $\\{0.001,0.01,0.05\\}$

    * –

Dropout rate: $\\{0,0.5\\}$

    * –

$K$: $10$

    * –

$\alpha$: $\\{0,0.1,0.2,0.5,0.7,0.9,1\\}$

    * –

dprate: $\\{0,0.5,0.7\\}$

##  Appendix F Additional Experiments

###  F.1 Experimental Results on Homophilic Benchmark Datasets

Competitive Performance on Real-World Homophilic Datasets. Table 3 summarizes
the averaged accuracy (the micro-F1 score) and standard deviation of semi-
supervised node classification on homophilic benchmark datasets. Table 3 shows
that the performance of pGNN is very close to APPNP, JKNet, GCN on Cora,
CiteSeer, PubMed datasets and slightly outperforms all baselines on Computers,
Photo, CS, Physics datasets. Moreover, we observe that pGNNs outperform GPRGNN
on all homophilic datasets, which confirms that pGNNs work better under weak
supervised information (2.5% training rate) as discussed in Remark 3. We also
see that all GNN models work significantly better than MLP on all homophilic
datasets. It illustrates that the graph topological information is helpful for
the label prediction tasks. Notably, 1.0GNN is slightly worse than the other
pGNNs with larger $p$, which suggests to use $p\approx 2$ for homophilic
graphs. Overall, the results of Table 3 indicates that pGNNs obtain
competitive performance against all baselines on homophilic datasets.

Table 3: Results on homophilic benchmark datasets. Averaged accuracy (%) for 100 runs. Best results are outlined in bold and the results within $95\%$ confidence interval of the best results are outlined in underlined bold. OOM denotes out of memory. Method | Cora | CiteSeer | PubMed | Computers | Photo | CS | Physics
---|---|---|---|---|---|---|---
MLP | $43.47_{\pm 3.82}$ | $46.95_{\pm 2.15}$ | $78.95_{\pm 0.49}$ | $66.11_{\pm 2.70}$ | $76.44_{\pm 2.83}$ | $86.24_{\pm 1.43}$ | $92.58_{\pm 0.83}$
GCN | $76.23_{\pm 0.79}$ | $62.43_{\pm 0.81}$ | $83.72_{\pm 0.27}$ | $84.17_{\pm 0.59}$ | $90.46_{\pm 0.48}$ | $90.33_{\pm 0.36}$ | $94.46_{\pm 0.08}$
SGC | $77.19_{\pm 1.47}$ | $\bm{64.10}_{\pm 1.36}$ | $79.26_{\pm 0.69}$ | $84.32_{\pm 0.59}$ | $89.81_{\pm 0.57}$ | $91.06_{\pm 0.05}$ | OOM
GAT | $75.62_{\pm 1.01}$ | $61.28_{\pm 1.09}$ | $83.60_{\pm 0.22}$ | $82.72_{\pm 1.29}$ | $90.48_{\pm 0.57}$ | $89.96_{\pm 0.27}$ | $93.96_{\pm 0.21}$
JKNet | $77.19_{\pm 0.98}$ | $63.32_{\pm 0.95}$ | $82.54_{\pm 0.43}$ | $79.94_{\pm 2.47}$ | $88.29_{\pm 1.64}$ | $89.69_{\pm 0.66}$ | $93.92_{\pm 0.32}$
APPNP | $\bm{79.58}_{\pm 0.59}$ | $63.02_{\pm 1.10}$ | $\bm{84.80}_{\pm 0.22}$ | $83.32_{\pm 1.11}$ | $90.42_{\pm 0.53}$ | $91.54_{\pm 0.24}$ | $\bm{94.93}_{\pm 0.06}$
GPRGNN | $76.10_{\pm 1.30}$ | $61.60_{\pm 1.69}$ | $83.16_{\pm 0.84}$ | $82.78_{\pm 1.87}$ | $89.81_{\pm 0.66}$ | $90.59_{\pm 0.38}$ | $\bm{\underline{94.72}}_{\pm 0.16}$
1.0GNN | $77.59_{\pm 0.69}$ | $63.19_{\pm 0.98}$ | $83.21_{\pm 0.30}$ | $84.46_{\pm 0.89}$ | $\bm{\underline{90.69}}_{\pm 0.66}$ | $91.46_{\pm 0.50}$ | $\bm{\underline{94.72}}_{\pm 0.37}$
1.5GNN | $78.86_{\pm 0.75}$ | $\bm{\underline{63.80}}_{\pm 0.79}$ | $83.65_{\pm 0.17}$ | $\bm{85.03}_{\pm 0.90}$ | $\bm{90.91}_{\pm 0.50}$ | $\bm{\underline{92.12}}_{\pm 0.40}$ | $\bm{\underline{94.90}}_{\pm 0.16}$
2.0GNN | $78.93_{\pm 0.60}$ | $\bm{\underline{63.65}}_{\pm 1.08}$ | $84.19_{\pm 0.22}$ | $84.39_{\pm 0.85}$ | $90.40_{\pm 0.63}$ | $\bm{92.28}_{\pm 0.47}$ | $\bm{94.93}_{\pm 0.14}$
2.5GNN | $78.87_{\pm 0.57}$ | $63.28_{\pm 0.97}$ | $\bm{\underline{84.45}}_{\pm 0.18}$ | $83.85_{\pm 0.87}$ | $89.82_{\pm 0.64}$ | $91.94_{\pm 0.40}$ | $\bm{\underline{94.87}}_{\pm 0.11}$

###  F.2 Experimental Results of Aggregation Weight Entropy Distribution

Here we present the visualization results of the learned aggregation weight
entropy distribution of pGNNs and GAT on all benchmark datasets. Fig. 5 and
Fig. 6 show the results obtained on homophilic and heterophilic benchmark
datasets, respectively.

Figure 5: Aggregation weight entropy distribution of homophilic benchmark
graphs. Low entropy means high degree of concentration, vice versa. An entropy
of zero means all aggregation weights are on one source node.

Figure 6: Aggregation weight entropy distribution of heterophilic benchmark
graphs. Low entropy means high degree of concentration and vice versa. An
entropy of zero means all aggregation weights are on one source node.

We observe from Fig. 5 that the aggregation weight entropy distributions
learned by pGNNs and GAT on homophilic benchmark datasets are similar to the
uniform cases, which indicates that aggregating and transforming node features
over the original graph topology is very helpful for label prediction. It
explains why pGNNs and GNN baselines obtained similar performance on
homophilic benchmark datasets and all GNN models significantly outperform MLP.

Contradict to the results on homophilic graphs shown in Fig. 5, Fig. 6 shows
that the aggregation weight entropy distributions of pGNNs on heterophilic
benchmark datasets are very different from that of GAT and the uniform cases.
We observe from Fig. 6 that the entropy of most of the aggregation weights
learned by pGNNs are around zero, which means that most aggregation weights
are on one source node. It indicates that the graph topological information in
these heterophilic benchmark graphs is not helpful for label prediction.
Therefore, propagating and transforming node features over the graph topology
could lead to worse performance than MLPs, which validates the results in
Table 3 that the performance of MLP is significantly better most GNN baselines
on all heterophilic graphs and closed to pGNNs.

###  F.3 Experimental Results on cSBM

In this section we present the experimental results on cSBM using sparse
splitting and dense splitting, respectively. We used the same settings in
Chien et al. (2021) in which the number of nodes $n=5000$, the number of
features $f=2000$, $\epsilon=3.25$ for all experiments. Table 4 reports the
results on cSBM with sparse splitting setting, which also are presented in
Fig. 2 and discussed in Sec. 5. Table 5 reports the results on cSBM with dense
splitting settings.

Table 4: Results on cSBM with sparse splitting setting. Average accuracy (%) for 20 runs. Best results are outlined in bold and the results within $95\%$ confidence interval of the best results are outlined in underlined bold. Method | $\phi=-1$ | $\phi=-0.75$ | $\phi=-0.5$ | $\phi=-0.25$ | $\phi=0$ | $\phi=0.25$ | $\phi=0.5$ | $\phi=0.75$ | $\phi=1$
---|---|---|---|---|---|---|---|---|---
MLP | $49.72_{\pm 0.36}$ | $51.42_{\pm 1.83}$ | $59.21_{\pm 1.01}$ | $\bm{61.57}_{\pm 0.38}$ | $\bm{\underline{61.70}}_{\pm 0.30}$ | $59.92_{\pm 1.88}$ | $57.20_{\pm 0.62}$ | $54.48_{\pm 0.48}$ | $50.09_{\pm 0.51}$
GCN | $57.24_{\pm 1.15}$ | $58.19_{\pm 1.46}$ | $57.30_{\pm 1.30}$ | $51.97_{\pm 0.44}$ | $54.45_{\pm 1.38}$ | $64.70_{\pm 2.38}$ | $82.45_{\pm 1.35}$ | $91.31_{\pm 0.54}$ | $76.07_{\pm 3.30}$
SGC | $55.98_{\pm 1.48}$ | $58.56_{\pm 1.40}$ | $56.97_{\pm 0.54}$ | $51.54_{\pm 0.22}$ | $52.69_{\pm 2.36}$ | $64.14_{\pm 1.05}$ | $79.88_{\pm 1.57}$ | $90.37_{\pm 0.09}$ | $75.94_{\pm 0.92}$
GAT | $59.72_{\pm 2.23}$ | $60.20_{\pm 2.14}$ | $55.38_{\pm 1.96}$ | $50.15_{\pm 0.55}$ | $53.05_{\pm 1.40}$ | $64.00_{\pm 2.03}$ | $81.04_{\pm 1.71}$ | $90.37_{\pm 1.33}$ | $78.24_{\pm 1.95}$
JKNet | $49.70_{\pm 0.39}$ | $49.75_{\pm 0.79}$ | $49.65_{\pm 0.52}$ | $48.93_{\pm 0.48}$ | $52.36_{\pm 2.09}$ | $62.76_{\pm 2.54}$ | $87.10_{\pm 1.52}$ | $\bm{\underline{97.43}}_{\pm 0.36}$ | $\bm{\underline{97.69}}_{\pm 0.52}$
APPNP | $48.45_{\pm 0.98}$ | $49.65_{\pm 0.46}$ | $53.31_{\pm 0.89}$ | $56.58_{\pm 0.58}$ | $60.10_{\pm 0.65}$ | $65.02_{\pm 2.23}$ | $82.95_{\pm 1.38}$ | $95.49_{\pm 0.43}$ | $89.85_{\pm 0.60}$
GPRGNN | $97.26_{\pm 0.66}$ | $94.81_{\pm 0.91}$ | $82.14_{\pm 0.47}$ | $61.15_{\pm 2.55}$ | $60.20_{\pm 0.76}$ | $62.90_{\pm 2.22}$ | $83.61_{\pm 1.28}$ | $96.96_{\pm 0.41}$ | $\bm{98.01}_{\pm 0.71}$
1.0GNN | $95.75_{\pm 1.21}$ | $93.06_{\pm 1.13}$ | $77.39_{\pm 4.21}$ | $\bm{\underline{61.38}}_{\pm 0.39}$ | $\bm{61.80}_{\pm 0.29}$ | $\bm{\underline{65.73}}_{\pm 2.11}$ | $85.85_{\pm 3.24}$ | $96.80_{\pm 0.87}$ | $97.40_{\pm 1.10}$