mathx"17

# Distributional Signals for Node Classification in Graph Neural Networks

Feng Ji, See Hian Lee, Kai Zhao, Wee Peng Tay, and Jielong Yang F. Ji, S. H.
Lee, K. Zhao and W. P. Tay are with the School of Electrical and Electronic
Engineering, Nanyang Technological University, 639798, Singapore. J. Yang is
with the School of Internet of Things Engineering, Jiangnan University, China.

###### Abstract

In graph neural networks (GNNs), both node features and labels are examples of
graph signals, a key notion in graph signal processing (GSP). While it is
common in GSP to impose signal smoothness constraints in learning and
estimation tasks, it is unclear how this can be done for discrete node labels.
We bridge this gap by introducing the concept of distributional graph signals.
In our framework, we work with the distributions of node labels instead of
their values and propose notions of smoothness and non-uniformity of such
distributional graph signals. We then propose a general regularization method
for GNNs that allows us to encode distributional smoothness and non-uniformity
of the model output in semi-supervised node classification tasks. Numerical
experiments demonstrate that our method can significantly improve the
performance of most base GNN models in different problem settings.

###### Index Terms:

Graph neural networks, graph signal processing, distributional signals, node
classification, smoothness

##  I Introduction

We consider the semi-supervised node classification problem [1] that
determines class labels of nodes in graphs given sample observations and
possibly node features. Numerous graph neural network (GNN) models have been
proposed to tackle this problem. One of the first models is the graph
convolutional network (GCN) [2]. Interpreted geometrically, a GCN aggregates
information such as node features from the neighborhood of each node of the
graph. Algebraically, this process is equivalent to applying a graph
convolution filter to node feature vectors. Subsequently, many GNN models with
different considerations are introduced. Popular models include the graph
attention network (GAT) [3] that learns weights between pairs of nodes during
aggregation, and the hyperbolic graph convolutional neural network (HGCN) [4]
that considers embedding of nodes of a graph in a hyperbolic space instead of
a Euclidean space. For inductive learning, GraphSAGE [5] is proposed to
generate low-dimensional vector representations for nodes that are useful for
graphs with rich node attribute information. While new models draw inspiration
from GCN, GCN itself is built upon the foundation of graph signal processing
(GSP).

GSP is a signal processing framework that handles graph-structured data [6, 7,
8]. A graph signal is a vector with each component corresponding to a node of
a graph. Examples include node features and node labels. Moreover,
convolutions used in models such as GCN are special cases of convolution
filters in GSP [6]. All these show the close connections between GSP theory
and GNNs.

In GSP, signal smoothness (over the graph) is widely used to regularize
inference tasks. Intuitively, a signal is smooth if its values are similar at
each pair of nodes connected by an edge. One popular way to formally define
signal smoothness is to use the Laplacian quadratic form. There are numerous
GSP tools that leverage a smooth prior of the graph signals. For example,
Laplacian (Tikhonov) regularization is proposed for noise removal in [6] and
signal interpolation [9]. In [10], it is used in graph signal in-painting and
anomaly detection. In [11], the same technique is used for graph topology
inference.

However, for GNNs, it is remarked in [12, Section 4.1.2] that “graph Laplacian
regularization can hardly provide extra information that existing GNNs cannot
capture”. Therefore, a regularization scheme based on feature propagation is
proposed. It is demonstrated to be effective by comparing with other methods
such as [13] and [14] based on adversarial learning and [15] that co-trains
GNN models with an additional agreement model, which gives the probability
that two nodes have the same label. We partially agree with the above
assertion regarding graph Laplacian regularization, while remaining
reservative about its full correctness. In this paper, we propose a method
that is inspired by Laplacian regularization. As our main contribution, we
introduce the notion of distributional graph signals, instead of considering
graph signals. Analogous to the graph signal smoothness defined using graph
Laplacian, we define the smoothness of distributional graph signals. Together
with another property known as non-uniformity, we devise a regularization
scheme for GNNs in node classification tasks. This approach is easy to
implement and can be used as a plug-in regularization term together with any
given base GNN model. Its effectiveness is demonstrated with numerical
results.

The rest of the paper is organized as follows. We formally introduce the
notion of distributional graph signals in Section II. In Section II-B, we
motivate by analyzing step graph signals naturally occurring in node
classification problems. We study two key properties, namely smoothness and
non-uniformity, of distributional graph signals in Section III-A and Section
III-B. Based on the discussions, we describe the proposed regularization
scheme in Section III-C. In Section III-D, we analyze the model and discuss
its relations with related works. We present experimental results in Section
IV and conclude in Section V.

##  II Distributional graph signals

In this section, we motivate and introduce distributional graph signals based
on GSP theory.

###  II-A GSP preliminaries and signal smoothness

In this subsection, we give a brief overview of GSP theory [6]. The focus is
on the discussion of graph signal smoothness.

Let ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ be an undirected graph with
${\mathcal{V}}$ the vertex set and ${\mathcal{E}}$ the edge set. Suppose the
size of the graph is $n=|{\mathcal{V}}|$. Fix an ordering of ${\mathcal{V}}$.
Then, the _space of graph signals_ can be identified with the vector space
$\mathbb{R}^{n}$, with a graph signal ${\bm{x}}\in\mathbb{R}^{n}$, which
assigns its $i$-th component to the $i$-th vertex of ${\mathcal{G}}$. By
convention, signals are in column form, and ${\bm{x}}(i)$ is the $i$-th
component of ${\bm{x}}$.

In GSP, the key notion is _the graph shift operator_. Though there are several
choices for the graph shift operator, in our paper, we consider a common
choice: ${\bm{L}}_{{\mathcal{G}}}$, the Laplacian of ${\mathcal{G}}$, defined
by
${\bm{L}}_{{\mathcal{G}}}={\bm{D}}_{{\mathcal{G}}}-{\bm{A}}_{{\mathcal{G}}}$,
where ${\bm{D}}_{{\mathcal{G}}},{\bm{A}}_{{\mathcal{G}}}$ are the degree
matrix and adjacency matrix of ${\mathcal{G}}$, respectively. The Laplacian
${\bm{L}}_{{\mathcal{G}}}$ is positive semi-definite and symmetric. By the
spectral theorem, it has an eigendecomposition
${\bm{L}}_{{\mathcal{G}}}={\bm{U}}_{{\mathcal{G}}}\Lambda_{{\mathcal{G}}}{\bm{U}}_{{\mathcal{G}}}^{\top}$.
In the decomposition, $\Lambda_{{\mathcal{G}}}$ is a diagonal matrix, whose
diagonal entries $\\{\lambda_{1},\ldots,\lambda_{n}\\}$ are eigenvalues of
${\bm{L}}_{{\mathcal{G}}}$. They are non-negative and we assume
$\lambda_{1}\leq\ldots\leq\lambda_{n}$. The associated eigenbasis
$\\{{\bm{u}}_{1},\ldots,{\bm{u}}_{n}\\}$ are the columns of
${\bm{U}}_{{\mathcal{G}}}$. In GSP, an eigenvector with a small eigenvalue
(and hence a small index) is considered to be smooth. The signal values of
such a vector have small fluctuations across the edges of ${\mathcal{G}}$.

Given a graph signal ${\bm{x}}$, its _graph Fourier transform_ is
$\hat{{\bm{x}}}={\bm{U}}_{{\mathcal{G}}}^{\top}{\bm{x}}$, or equivalently,
$\hat{{\bm{x}}}(i)=\langle{\bm{x}},{\bm{u}}_{i}\rangle$, for $1\leq i\leq n$.
The components $\hat{{\bm{x}}}(i)$ of $\hat{{\bm{x}}}$ are called the
_frequency components_ of ${\bm{x}}$. Same as above, the signal ${\bm{x}}$ is
smooth if $\hat{{\bm{x}}}(i)$ has a small absolute value for large $i$.
Quantitatively, we can define its _total variation_ by

$\displaystyle\mathcal{T}_{G}({\bm{x}})=\sum_{(v_{i},v_{j})\in{\mathcal{E}}}({\bm{x}}(i)-{\bm{x}}(j))^{2}={\bm{x}}^{\top}{\bm{L}}_{{\mathcal{G}}}{\bm{x}}.$
(1)

It is straightforward to compute that
$\mathcal{T}_{G}({\bm{u}}_{i})=\lambda_{i}$. This observation indicates that
it is reasonable to use total variation as a measure of smoothness. Minimizing
the total variation of graph signals has many applications in GSP as we have
pointed out in Section I.

###  II-B Step graph signals

Let ${\mathbb{S}}$ be a finite set of numbers. A _step graph signal with
respect to (w.r.t.) ${\mathbb{S}}$_ is a graph signal ${\bm{x}}$ such that all
its components take values in ${\mathbb{S}}$, i.e.,
${\bm{x}}\in{\mathbb{S}}^{n}$.

######  Example 1.

For the simplest example, consider the classical Heaviside function $H$ on
$\mathbb{R}$ defined by $H(x)=1$, for $x>0$ and $H(x)=0$, for $x\leq 0$. It is
a non-smooth function, as it is not even continuous at $x=0$. On the other
hand, let ${\mathcal{G}}$ be the path graph with $2m+1$ nodes embedded on the
real line by identifying the nodes of ${\mathcal{G}}$ with the integers in the
interval $[-m,m]$. Then $H$ induces a step graph signal ${\bm{h}}$ on
${\mathcal{G}}$. Same as the Heaviside function $H$, the signal ${\bm{h}}$
should be considered to be a non-smooth graph signal.

Step graph signals occur naturally in semi-supervised node classification
tasks. In particular, if ${\mathbb{S}}$ is the set of all possible class
labels, then the labels of the nodes of ${\mathcal{G}}$ form a step graph
signal ${\bm{c}}$ w.r.t. ${\mathbb{S}}$ on ${\mathcal{G}}$. We expect that
analogous to Example 1, ${\bm{c}}$ can possibly be non-smooth. To demonstrate,
we analyze ${\bm{c}}$ using its Fourier transform $\hat{{\bm{c}}}$ (cf.
Section II-A). More specifically, we take ${\mathcal{G}}_{0}$ to be the main
connected component of the Cora graph [16] and ${\bm{c}}$ to be the ground
truth labels. We also generate a random signal ${\bm{r}}$ following the same
empirical distribution (on ${\mathbb{S}}$) estimated using ${\bm{c}}$. We show
plots of both Fourier transforms $\hat{{\bm{c}}}$ and $\hat{{\bm{r}}}$ in Fig.
1. We see that the high-frequency components of the ground truth labels
${\bm{c}}$ can be large, and it is even possible that its spectrum, i.e.,
frequency components, resemble that of a random signal. The observations
support our speculation about the non-smoothness of the step signals.
Therefore, in order to leverage signal smoothness to enhance model
performance, we need to find an alternative to the step (label) signals. This
also supports the remark of [12] regarding Laplacian regularization from a
different point of view (see also the experiments in Section IV-B).

Figure 1: Plots of $\hat{{\bm{c}}}$ and $\hat{{\bm{r}}}$ for Cora.

###  II-C Distributional graph signals and marginals

We want to use probability theory to introduce the notion of “smoothness” for
step signals in the next section. For preparation, in this subsection, we
formally introduce distributional graph signals and discuss how they arise
naturally in GNNs. Our theory (e.g., in Definition 1 and equation 3) is based
on probability measures defined on metric spaces. To this end, we endow any
discrete space ${\mathbb{S}}$ with the discrete metric $d(s_{1},s_{2})=1$ if
$s_{1}\neq s_{2}\in{\mathbb{S}}$ and $0$ otherwise. For $\mathbb{R}^{n}$, one
can use the usual Euclidean metric or any other norm.

######  Definition 1.

For a metric space ${\mathbb{M}}={\mathbb{S}}^{n}$ or $\mathbb{R}^{n}$, let
$\mathcal{P}({\mathbb{M}})$ be the space of probability measures on
${\mathbb{M}}$ (w.r.t. the Borel $\sigma$-algebra) having finite second
moments. An element $\mu\in\mathcal{P}({\mathbb{M}})$ is called a
_distributional graph signal_. The marginals of a distributional graph signal
$\mu$ are the marginal distributions $\mathscr{N}=\\{\mu_{i},1\leq i\leq n\\}$
w.r.t. the $n$ coordinates of either $\mathbb{S}^{n}$ or $\mathbb{R}^{n}$.

To understand why a distribution $\mu$ is related to GSP, consider a step
graph signal (resp. ordinary graph signal) ${\bm{x}}$ in ${\mathbb{S}}^{n}$
(resp. $\mathbb{R}^{n}$). It induces the delta distribution
$\delta_{{\bm{x}}}$ in $\mathcal{P}({\mathbb{S}}^{n})$ supported at
${\bm{x}}$. Its marginals are the delta distributions
$\\{\delta_{{\bm{x}}(i)},1\leq i\leq n\\}$. Therefore, Definition 1 subsumes
ordinary graph signals as special cases.

Figure 2: We give an example of the marginals of a distributional graph
signal. We notice that for ${\bm{X}}_{1,:}$, the probability weights are equal
and it is hard to determine the $\operatorname*{arg\,max}$ of the $3$ classes.
Moreover, the weights of the $2$nd class have large differences along a few
edges, e.g., $(v_{2},v_{3}),(v_{3},v_{4})$. These features make prediction
unreliable. We propose solutions to such issues in this paper.

Distributional graph signals also occur naturally in GNN models (illustrated
in Fig. 2). As in the previous subsection, let ${\mathbb{S}}$ be the set of
all possible class labels of size $m$ and $M_{n,m}(\mathbb{R})$ be the space
of $n\times m$ matrices. For a typical model $\mathcal{M}$ such as the GCN,
the last stage of the pipeline usually consists of the following steps111To
simplify the presentation, we assume that all models considered have the
intermediate softmax step, though some implementations omit this part using
the fact that the exponential function is increasing.:

$\displaystyle\begin{split}&\text{Logits: }{\bm{O}}\in
M_{n,m}(\mathbb{R})\xrightarrow{\text{Softmax}:\phi}{\bm{X}}=\phi({\bm{O}})\in
M_{n,m}(\mathbb{R})\\\ &\xrightarrow{\operatorname*{arg\,max}}\text{Output
labels: }{\bm{c}}\in{\mathbb{S}}^{n}.\end{split}$ (2)

The $i$-th row ${\bm{X}}_{i,:}$ of ${\bm{X}}$ can be viewed as the weights of
a probability distribution $\mu_{{\bm{X}},i}$ on ${\mathbb{S}}$. They do not
directly give a distributional graph signal in
$\mathcal{P}({\mathbb{S}}^{n})$. However, we shall interpret
$\mathscr{N}_{{\bm{X}}}=\\{\mu_{{\bm{X}},i},1\leq i\leq n\\}$ as the marginals
of some unknown distributional graph signal, and $\mathscr{N}_{{\bm{X}}}$ is
the main subject of study in this paper. In order to mimic the ways smoothness
of graph signals is used in GSP, we introduce an appropriate notion of the
total variation of distributional graph signals in the next section. Moreover,
we also want to study a unique property of the non-uniformity of
distributional graph signals. A signal processing theory of distributional
graph signal is developed in [17].

##  III Regularized graph neural networks

In this section, we study smoothness and non-uniformity properties of
distributional graph signals in Section III-A and Section III-B respectively.
Smoothness refers to the similarity of distributions for neighboring nodes,
which is different from the classical point of view. Moreover, it only makes
sense to talk about non-uniformity when dealing with distributional graph
signals. Each of these subsections yields an expression that we want to
minimize. They are combined in Section III-C to give the proposed
regularization term. Both smoothness and non-uniformity are important for the
effectiveness of the proposed model as we shall demonstrate with numerical
experiments.

###  III-A Smoothness of distributional graph signals

The goal of this subsection is to introduce and compare different notions of
total variation associated with distributional graph signals and their
marginals that leverage signal smoothness. First of all, given
$\mu\in\mathcal{P}({\mathbb{M}})$ for ${\mathbb{M}}={\mathbb{S}}^{n}$ or
$\mathbb{R}^{n}$, its total variation can be modified directly from (1) as
follows:

$\displaystyle\mathcal{T}_{G}(\mu)=\mathbb{E}_{{\bm{x}}\sim\mu}\mathcal{T}_{G}({\bm{x}})=\int\sum_{(v_{i},v_{j})\in{\mathcal{E}}}d({\bm{x}}(i),{\bm{x}}(j))^{2}\,\mathrm{d}\mu({\bm{x}}).$
(3)

where $d$ is the metric on ${\mathbb{S}}$ or $\mathbb{R}$.

However, in many cases of interest such as GNNs, only marginals of some $\mu$
are observed (cf. Section II-C). We want to define total variation in such a
situation as well. We borrow ideas from the prototype of the Wasserstein
metric [18], which we recall now. The notations and assumptions follow those
of Definition 1.

######  Definition 2.

For $\mu_{1},\mu_{2}\in\mathcal{P}({\mathbb{M}})$, the _Wasserstein metric_
$W(\mu_{1},\mu_{2})$ between $\mu_{1},\mu_{2}$ is defined by

$\displaystyle
W(\mu_{1},\mu_{2})^{2}=\inf_{\gamma\in\Gamma(\mu_{1},\mu_{2})}\int
d(x,y)^{2}\,\mathrm{d}\gamma(x,y),$

where $\Gamma(\mu_{1},\mu_{2})$ is the set of _couplings_ of
$\mu_{1},\mu_{2}$, i.e., the collection of probability measures on
${\mathbb{M}}\times{\mathbb{M}}$ whose marginals are $\mu_{1}$ and $\mu_{2}$,
respectively.

It can be verified that $W(\cdot,\cdot)$ indeed defines a metric on
$\mathcal{P}({\mathbb{M}})$ (see [18]). As a special case, if $\delta_{x}$ and
$\delta_{y}$ are delta distributions supported on $x,y\in{\mathbb{M}}$, then
$W(\delta_{x},\delta_{y})=d(x,y)$. The key insight is that we want to take
infimum over all possible distributions given the prescribed marginals; and
whatever we define, it should subsume (1) as a special case for a collection
of delta distributions.

######  Definition 3.

Given $\mathscr{N}=\\{\mu_{i},1\leq i\leq n\\}$ with
$\mu_{i}\in\mathcal{P}({\mathbb{S}})$ (resp.
$\mu_{i}\in\mathcal{P}(\mathbb{R})$), for $1\leq i\leq n$, then the total
variation of $\mathscr{N}$ is defined as:

$\displaystyle\mathcal{T}_{G}(\mathscr{N})=\inf_{\mu\in\Gamma(\mathscr{N})}\mathcal{T}_{G}(\mu),$
(4)

where $\Gamma(\mathscr{N})$ is the collection of all distributional graph
signals in $\mathcal{P}({\mathbb{S}}^{n})$ (resp.
$\mathcal{P}(\mathbb{R}^{n})$) whose marginals agree with $\mathscr{N}$.

Though an important theoretical tool, the Wasserstein metric is usually
difficult to compute explicitly. On the other hand, if ${\mathcal{G}}$ is the
graph with $2$ nodes connected by an edge and
$\mathscr{N}=\\{\mu_{1},\mu_{2}\\}$, then
$\mathcal{T}_{G}(\mathscr{N})=W(\mu_{1},\mu_{2})^{2}$. As a consequence,
finding the exact value of $\mathcal{T}_{G}(\mathscr{N})$ can be challenging.
Therefore, we next introduce approximations that can be more readily
estimated.

In graph signal processing, a useful strategy to understand a general graph is
to study its spanning trees [19, 20]. A set
${\mathbb{T}}=\\{T_{1},\ldots,T_{k}\\}$ of $k$-spanning trees of $G$ _covers_
if each edge of $G$ is contained in $T_{i}$ for some $1\leq i\leq k$.

######  Definition 4.

Define the _covering total variation_ $\mathcal{T}^{c}_{G}(\mathscr{N})$ by

$\displaystyle\mathcal{T}^{c}_{G}(\mathscr{N})=\inf_{{\mathbb{T}}=\\{T_{1},\ldots,T_{k}\\}}\sum_{1\leq
i\leq k}\mathcal{T}_{T_{i}}(\mathscr{N}),$

where the infimum is over all covers $\mathbb{T}$ of $G$ by spanning trees.

$\mathcal{T}_{G}(\mathscr{N})$ is the most principled notion of total
variation. The proxy $\mathcal{T}^{c}_{G}(\mathscr{N})$ retains important
features of $\mathcal{T}_{G}(\mathscr{N})$ that encode edgewise differences of
distributional signals while taking the entire graph $G$ into consideration.
Moreover, $\mathcal{T}^{c}_{G}(\mathscr{N})$ can be readily estimated using
more familiar quantities.

For $\mathscr{N}=\\{\mu_{i},1\leq i\leq n\\}$ with
$\mu_{i}\in\mathcal{P}({\mathbb{S}})$, let $\bm{\mu}_{s}\in\mathbb{R}^{n}$ be
the graph signal of probability weights of $s\in{\mathbb{S}}$, i.e.,
$\bm{\mu}_{s}(i)$ is the probability weight of $\mu_{i}$ at $s$. We can also
stack these signals as a matrix ${\bm{X}}_{\mathscr{N}}$ whose columns are
$\bm{\mu}_{s}$. Based on the GSP version of total variation (cf. 1), we
introduce two more versions of total variation that are easy to compute.

######  Definition 5.

Given $\mathscr{N}=\\{\mu_{i},1\leq i\leq n\\}$ with
$\mu_{i}\in\mathcal{P}({\mathbb{S}})$, we define the $\ell^{1}$ and $\ell^{2}$
versions of total variations as:

  * •

$\mathcal{T}_{G,1}(\mathscr{N})=\sum_{s\in{\mathbb{S}}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}|\bm{\mu}_{s}(i)-\bm{\mu}_{s}(j)|$,
and

  * •

$\mathcal{T}_{G,2}(\mathscr{N})=\sum_{s\in{\mathbb{S}}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}\big{(}\bm{\mu}_{s}(i)-\bm{\mu}_{s}(j)\big{)}^{2}=\operatorname{Tr}({\bm{X}}_{\mathscr{N}}^{\top}{\bm{L}}_{{\mathcal{G}}}{\bm{X}}_{\mathscr{N}})$,
where $\operatorname{Tr}$ is the matrix trace.

The complexity of computing either $\mathcal{T}_{G,1}$ or $\mathcal{T}_{G,2}$
is at most $\mathcal{O}(|{\mathbb{S}}||{\mathcal{E}}|)$, and it involves only
matrix multiplication for $\mathcal{T}_{G,2}$. In addition,
$\mathcal{T}^{c}_{G},\mathcal{T}_{G,1}$ and $\mathcal{T}_{G,2}$ satisfy the
following relation.

######  Theorem 1.

Given $\mathscr{N}=\\{\mu_{i},1\leq i\leq n\\}$ with
$\mu_{i}\in\mathcal{P}({\mathbb{S}})$, we have

$\displaystyle\mathcal{T}_{G,2}(\mathscr{N})\leq\mathcal{T}_{G,1}(\mathscr{N})$
$\displaystyle\leq
2\min\big{(}\mathcal{T}^{c}_{G}(\mathscr{N}),\mathcal{T}_{G}(\mathscr{N})\big{)}$
$\displaystyle\leq c_{1}\mathcal{T}_{G,1}(\mathscr{N})\leq
c_{2}\mathcal{T}_{G,2}(\mathscr{N})^{\frac{1}{2}},$

where $c_{1},c_{2}$ are constants independent of $\mathscr{N}$.

The discussion and proof of a more general result can be found in Appendix C.
The upshot is that if we expect $\mathcal{T}^{c}_{G}(\mathscr{N})$ to be small
for some $\mathscr{N}$, then so are necessarily
$\mathcal{T}_{G,1}(\mathscr{N})$ and $\mathcal{T}_{G,2}(\mathscr{N})$. In view
of computation cost, we mainly use $\mathcal{T}_{G,2}$ in the design of the
regularization model in Section III-C. It is discussed in many recent works
regarding unfavorable side effects of signal smoothness (cf. Section III-D).
However, our point of view is that the smoothness of “distributional graph
signals” might be a worthy alternative. This is conceptually different from
classical Laplacian-based regularization.

###  III-B Non-uniformity of distributional graph signals

As discussed in Section II-C, a base GNN model may output a matrix ${\bm{X}}$
(in (2)), with associated $\mathscr{N}_{{\bm{X}}}=\\{\mu_{{\bm{X}},i},1\leq
i\leq n\\}$. For node classification problems, it is desirable that there is
less ambiguity in the decision for each node so that one can pinpoint the
correct class label. Mathematically, this requires that each
$\mu_{{\bm{X}},i}$ deviates from the uniform distribution
$\mathcal{U}({\mathbb{S}})$ on the finite set of label classes ${\mathbb{S}}$,
measured by the Wasserstein metric
$W(\mu_{{\bm{X}},i},\mathcal{U}({\mathbb{S}}))$. The following result is
proved in Appendix C.

######  Lemma 1.

For a fixed sequence of non-positive numbers $(a_{i})_{1\leq i\leq n}$, there
is a constant $C$ independent of ${\bm{X}}$ such that

$\displaystyle\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})+C\geq
2\sum_{1\leq i\leq n}a_{i}W(\mu_{{\bm{X}},i},\mathcal{U}({\mathbb{S}}))^{2},$
(5)

where ${\bm{D}}$ is the diagonal matrix with diagonal entries $(a_{i})_{1\leq
i\leq n}$. Moreover,
$\operatorname{Tr}({\bm{X}}_{o}^{\top}{\bm{D}}{\bm{X}}_{o})\leq\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})\leq\operatorname{Tr}({\bm{X}}_{u}^{\top}{\bm{D}}{\bm{X}}_{u})$,
where ${\bm{X}}_{o}$ is a matrix with each row a one-hot vector and
${\bm{X}}_{u}$ is the matrix with each entry $1/|{\mathbb{S}}|$.

As we want each $\mu_{{\bm{X}},i}$ to deviate from the uniform distribution,
the right-hand side of 5 should be made small (as negative as possible). This
is ensured if the proxy $\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})$,
which is easy to compute, is small. Moreover, we notice that ${\bm{X}}_{u}$
(resp. ${\bm{X}}_{o}$) corresponds to uniform (resp. $\delta$) marginal
distributions. The second half of the statement suggests that minimizing
$\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})$ may drive marginals,
i.e., rows of ${\bm{X}}$, to comply with non-uniformity. Numerical experiments
in Section IV-B support that this proxy works well. We use this term in
conjunction with $\mathcal{T}_{G,2}$ introduced in Section III-A in our
proposed regularized model.

Regarding the numbers $(a_{i})_{1\leq i\leq n}$, we propose to assign a larger
number to a node with a larger degree. The reason is that we would like to
promote the non-uniformity of large degree nodes, as each such node can
influence more nodes during message passing.

###  III-C The loss function with regularization

Suppose we are given a base GNN model, denoted by $\mathcal{M}$. Let
${\bm{F}}$ be the matrix of input feature vectors and $\Theta$ be the
parameter space for $\mathcal{M}$. Assume $\mathcal{M}$ has a loss function
$\mathcal{L}_{\mathcal{M}}$, and the model is set to solve the optimization
problem $\min_{\theta\in\Theta}\mathcal{L}_{\mathcal{M}}({\bm{F}},\theta)$.

Consider the steps given in 2. We have ${\bm{X}}=\phi({\bm{O}})$ viewed as a
matrix of probability weights. For regularization, we introduce another loss
$\mathcal{L}_{0}$ to supplement $\mathcal{L}_{\mathcal{M}}$. Sections III-A
and III-B suggest that $\mathcal{L}_{0}$ should consist of two parts
$\mathcal{L}_{0}=\mathcal{L}_{1}+\mathcal{L}_{2}$. The loss $\mathcal{L}_{1}$
(cf. Definition 5) is related to the smoothness of the distributional graph
signal ${\bm{X}}$ and takes the form
$\mathcal{L}_{1}({\bm{X}})=\operatorname{Tr}({\bm{X}}^{\top}{\bm{L}}_{{\mathcal{G}}}{\bm{X}})=\operatorname{Tr}({\bm{X}}^{\top}({\bm{D}}_{{\mathcal{G}}}-{\bm{A}}_{{\mathcal{G}}}){\bm{X}})$.
On the other hand, $\mathcal{L}_{2}$ prevents the distributional signal from
being uniform and can be explicitly expressed as
$\mathcal{L}_{2}({\bm{X}})=\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})$
for a suitably chosen negative semi-definite diagonal matrix ${\bm{D}}$ (cf.
Lemma 1). Summing up $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$, we have
$\mathcal{L}_{0}({\bm{X}})=\operatorname{Tr}\big{(}{\bm{X}}^{\top}({\bm{D}}_{{\mathcal{G}}}-{\bm{A}}_{{\mathcal{G}}}+{\bm{D}}){\bm{X}}\big{)}$.
In our experiments, we take ${\bm{D}}={\bm{I}}_{n}-{\bm{D}}_{{\mathcal{G}}}$
and obtain the easily computable loss

$\displaystyle\mathcal{L}_{0}({\bm{X}})=\operatorname{Tr}\big{(}{\bm{X}}^{\top}({\bm{I}}_{n}-{\bm{A}}_{{\mathcal{G}}}){\bm{X}}\big{)}.$

In summary, if we express the output of the second last layer of $\mathcal{M}$
as ${\bm{O}}=\psi({\bm{F}},\theta)$ of input ${\bm{F}}$ and model parameters
$\theta$, then the regularized model R-$\mathcal{M}$ of $\mathcal{M}$ solves
the optimization:

$\displaystyle\min_{\theta\in\Theta}\mathcal{L}_{\text{R-}\mathcal{M}}({\bm{F}},\theta)=\min_{\theta\in\Theta}\mathcal{L}_{\mathcal{M}}({\bm{F}},\theta)+\eta\cdot\mathcal{L}_{0}\big{(}\phi\circ\psi({\bm{F}},\theta)\big{)},$
(6)

where the coefficient $\eta$ is a tunable hyperparameter and $\circ$ denotes
function composition. In the regularized model R-$\mathcal{M}$, we do not make
any other changes to the based model $\mathcal{M}$ apart from using the new
loss function $\mathcal{L}_{\text{R-}\mathcal{M}}$ during training. A
schematic illustration is shown in Fig. 3 (a).

Figure 3: This figure illustrates in (a) how R-$\mathcal{M}$ is constructed
upon the base model $\mathcal{M}$, and in (b) an R-module that can be used at
different places in GNN models.

###  III-D Further discussions

Most regularization methods introduce a penalty to supplement the base model
loss. For example, in the recent work, P-reg [12] proposes to apply a
propagation matrix (e.g., the normalized adjacency matrix) to the output
features. The model penalizes large discrepancies, measured by a metric such
as squared error distance, between the original and transformed features. In
LEreg [21], intra-energy and inter-energy losses are introduced. Both are
variants of the total variation 1. The novelty lies in introducing a “merged
graph” with each node representing a whole label class. All these works are
related to the smoothness prior used in GSP theory (Section II-A). In this
paper, we take a fundamentally different view of graph signals by _treating a
distribution as a signal_. This allows a _principled_ and _straightforward_
adaption of existing GSP approaches, as long as we have suitable notions of
the signal’s total variation. Moreover, non-uniformity does not have a
counterpart for ordinary graph signals, which are essentially delta
distributions.

Compared to GSP as well as other regularization methods such as P-reg and
LEreg, a salient feature of our method is that the matrix
${\bm{I}}_{n}-{\bm{A}}_{{\mathcal{G}}}$ in $\mathcal{L}_{0}$ is in general not
positive semi-definite. Therefore, the infimum of $\mathcal{L}_{0}({\bm{X}})$
can be $-\infty$ if the domain of the entries of ${\bm{X}}$ is unbounded. This
is another reason why restricting to distributional graph signals is
essential.

The proposed regularization is primarily for node classification as
distributional graph signals can be interpreted as the likelihoods of class
labels. However, the method can be extended to other tasks through an _R-
module_ (Fig. 3 (b)) that consists of the following steps:

  * •

Apply $\phi$ (e.g., $\mathrm{softmax}$) that turns a feature ${\bm{X}}$ into a
matrix of probability weights $\phi({\bm{X}})$.

  * •

Plug $\phi({\bm{X}})$ in the loss $\eta_{{\bm{X}}}\mathcal{L}_{0}$ with
tunable coefficient $\eta_{{\bm{X}}}$.

Given a base model, multiple R-modules can be inserted at different places of
the model pipeline, and all the losses are combined with the original loss of
the model. The insight is that useful node features for graph learning tasks
should be inherently associated with smooth and non-uniform distributional
graph signals. We demonstrate this approach with link prediction and graph
classification tasks in Appendix B.

_More on related works._ Though it is argued in [12] and supported by our
evidence in Section II-B that Laplacian regularization may have its drawbacks,
it has achieved a certain degree of success in earlier works [22, 23, 24].
These methods are based on the insight that neighboring nodes are likely to
have the same labels. In the language of GSP, they intend to leverage the
smoothness of the class label graph signal.

On the other hand, though graph signal smoothness has been fundamental in both
GSP and GNNs, the negative effect of over-smoothing has also been examined
[25, 26]) and models are proposed to alleviate it. For example, apart from the
more recent works such as [12, 21] that have already been discussed in detail,
PairNorm proposed in [27] encourages the similarity between connected nodes,
and at the same time adds a negative term based on distances between
disconnected pairs. MADReg in [28] proposes to use step size limits to make
the graph nodes receive less interference noise. In [29], randomly dropping
nodes is proposed to reduce the convergence speed of over-smoothing. Adding
skip connections is also introduced in [30, 31].

In our paper, we take a different point of view by not considering the
smoothness of ordinary graph signals. Instead, we speculate that properties
such as smoothness and non-uniformity of distributional graph signals may play
important roles. Moreover, requiring a distributional graph signal to satisfy
non-uniformity partially prevents the unfavorable situation that many
connected nodes have similar marginal distributions that are approximately
uniform. This can be viewed as a countermeasure to over-smoothing
intrinsically contained in our approach.

##  IV Experiments

In this section, we verify the empirical performance of our proposed
regularization method based on both Euclidean and hyperbolic GNN models. We
consider node classification under both transductive and inductive settings.
The datasets used include Cora, Citeseer, Pubmed, (Amazon) Photo, CS, Airport,
Disease, and PPI [16, 32, 33, 34, 35, 4, 36]. Their statistics are given in
Appendix A. No further tuning of the base model is needed. We compare with
base models and benchmarks in Section IV-A, and with variants of our model in
Section IV-B. Other graph learning tasks are discussed in Appendix B.

###  IV-A Performance on node classification problems

####  IV-A1 Transductive learning models

In this subsection, we consider transductive tasks, in which there is a single
graph containing both labeled training nodes and test nodes. The base models
include GCN [2], GAT [3], GraphSAGE [5] (abbreviated as SAGE) and GraphCON
[37] (abbreviated as CON).

In addition to comparison with base models, we use models having a similar
structure to our approach as benchmarks. More specifically, we implement
different versions of P-reg in [12], denoted by P-GCN and P-GAT (based on
popular GNN models GCN, GAT), and different versions of LEReg [21], denoted by
L-GCN and L-GAT. The parameters are tuned as suggested by the respective
papers. We also have the Laplacian method, denoted by LAP, in which the final
class label ${\bm{c}}$ is used in the loss $\mathcal{L}_{0}$. For a fair
comparison, tests are performed under the same hardware and software
environment. Similar considerations are applied in Section IV-A2 and IV-A3.

We also compare with BVAT [14], GAM [15], and GraphAT [13] (the results are
taken from literatures though they are unreported for Photo and CS). The test
accuracies (in $\%$) are shown in Table I. In each row, best performers are
highlighted in  blue and  red for our approach and benchmarks, respectively.
An underlined entry means no noticeable performance improvement over the base
model is observed. In general, we see that our proposed regularized models
improve upon their respective base models with significant performance gain in
many cases. Moreover, our method can match up with or even outperform many
benchmarks.

TABLE I: Transductive learning models.

| Comparison w/ base models | Benchmarks
---|---|---
Dataset | GCN | R-GCN | GAT | R-GAT | SAGE | R-SAGE | CON |  R-CON | LAP | P-GCN | P-GAT | L-GCN | L-GAT | BVAT | GAM | GraphAT
Cora | $81.0\atop\pm 1.07$ | $83.4\atop\pm 0.24$ | $83.1\atop\pm 0.61$ | $83.7\atop\pm 0.40$ | $81.8\atop\pm 0.65$ | $82.6\atop\pm 0.31$ | $83.9\atop\pm 1.12$ |  $\bm{84.4}\atop\pm 1.52$ | $80.7\atop\pm 0.89$ | $80.4\atop\pm 0.61$ | $81.2\atop\pm 0.65$ | $81.5\atop\pm 0.45$ | $82.6\atop\pm 0.75$ | $\bm{83.4}\atop\pm 1.01$ | $82.3\atop\pm 0.48$ | $82.5\atop\pm 0.60$
Citeseer | $70.6\atop\pm 0.48$ | $73.6\atop\pm 0.70$ | $68.7\atop\pm 0.69$ | $70.7\atop\pm 1.23$ | $70.6\atop\pm 0.47$ | $72.0\atop\pm 0.28$ | $73.5\atop\pm 1.27$ |  $\bm{74.6}\atop\pm 1.16$ | $70.7\atop\pm 0.64$ | $70.4\atop\pm 0.87$ | $70.3\atop\pm 0.76$ | $71.1\atop\pm 0.59$ | $69.0\atop\pm 0.94$ | $\bm{73.9}\atop\pm 0.54$ | $72.7\atop\pm 0.62$ | $73.5\atop\pm 0.38$
Pubmed | $79.2\atop\pm 0.29$ | $79.9\atop\pm 0.40$ | $76.9\atop\pm 0.77$ | $\underline{76.9}\atop\pm 0.77$ | $77.8\atop\pm 0.31$ | $78.7\atop\pm 0.53$ | $79.1\atop\pm 1.16$ |  $\bm{80.0}\atop\pm 1.18$ | $79.1\atop\pm 0.73$ | $79.0\atop\pm 0.31$ | $76.9\atop\pm 0.95$ | $78.7\atop\pm 0.25$ | $73.5\atop\pm 0.47$ | $78.0\atop\pm 0.84$ | $\bm{79.6}\atop\pm 0.63$ | $79.1\atop\pm 0.20$
Photo | $91.1\atop\pm 0.55$ | $91.4\atop\pm 0.49$ | $91.4\atop\pm 0.41$ | $\bm{92.5}\atop\pm 0.42$ | $90.0\atop\pm 0.36$ | $92.2\atop\pm 0.42$ | $90.3\atop\pm 0.64$ |  $91.1\atop\pm 0.43$ | $91.1\atop\pm 0.43$ | $91.6\atop\pm 0.20$ | $91.6\atop\pm 1.10$ | $91.5\atop\pm 0.28$ | $\bm{92.0}\atop\pm 0.45$ | - | - | -
CS | $90.5\atop\pm 0.14$ | $91.2\atop\pm 0.65$ | $90.8\atop\pm 0.40$ | $91.7\atop\pm 0.25$ | $90.5\atop\pm 0.08$ | $\bm{91.8}\atop\pm 0.10$ | $90.5\atop\pm 0.44$ |  $90.8\atop\pm 0.36$ | $90.6\atop\pm 0.11$ | $90.5\atop\pm 0.14$ | $90.8\atop\pm 0.40$ | $90.5\atop\pm 0.10$ | $\bm{90.8}\atop\pm 0.34$ | - | - | -

####  IV-A2 Hyperbolic models

Base models considered in Section IV-A1 generate embedding of nodes in
Euclidean spaces. However, it is argued in works such as [4] that for certain
datasets such as Airport and Disease, one should consider embedding nodes in
hyperbolic spaces and perform feature aggregation in the tangent spaces of
hyperbolic spaces. Such a consideration is plausible as certain graphs are
inherently hyperbolic (measured by $\delta$-hyperbolicity, see [38]).

In this subsection for Airport and Disease datasets, we use hyperbolic
versions of their Euclidean counterparts HGCN [4], and HGAT [39] as base
models. We also consider the interactive model GIL that combines both
Euclidean and hyperbolic approaches [40]. The comparison results are shown in
Table II. Again, we see a general improvement by using the proposed
regularization, which yields performance comparable with benchmarks.

TABLE II: Hyperbolic models.

| Comparison w/ base models | Benchmarks
---|---|---
Dataset | HGCN | R-HGCN | HGAT | R-HGAT | GIL |  R-GIL | LAP | P-HGCN | P-HGAT | L-HGCN | L-HGAT
Airport | $88.8\atop\pm 1.38$ | $89.3\atop\pm 1.42$ | $89.9\atop\pm 0.83$ | $91.0\atop\pm 1.63$ | $90.3\atop\pm 2.12$ |  $\bm{91.4}\atop\pm 1.39$ | $88.8\atop\pm 1.38$ | $88.9\atop\pm 1.42$ | $\bm{90.0}\atop\pm 1.11$ | $88.7\atop\pm 1.57$ | $89.7\atop\pm 1.14$
Disease | $91.3\atop\pm 1.04$ | $\bm{92.1}\atop\pm 1.09$ | $90.4\atop\pm 1.24$ | $91.3\atop\pm 1.31$ | $91.2\atop\pm 1.52$ |  $\underline{91.2}\atop\pm 1.52$ | $92.2\atop\pm 1.01$ | $92.1\atop\pm 1.07$ | $90.2\atop\pm 2.14$ | $\bm{92.2}\atop\pm 0.95$ | $90.4\atop\pm 1.53$

####  IV-A3 Inductive learning models

In contrast with transductive learning, inductive learning requires one to
deal with unseen data outside the training set. For example, in the PPI
dataset, different graphs correspond to different human tissues and we test
two out of 24 graphs that are unseen during training. Though the citation
datasets Cora, Citeseer, and Pubmed are for transductive learning, we modify
the datasets following [41] that use the induced subgraph of training nodes
during training. During validation and testing, we use induced subgraphs of
validation nodes and testing nodes with training nodes, respectively. The
nodes that are being predicted are unseen during training.

The base models are GCN, GAT, as well as GraphSAGE that is primarily designed
for inductive learning tasks. The performance comparison between the base
models and their regularized versions is shown in Table III. The conclusion
agrees with those observed in the previous subsections.

TABLE III: Inductive learning models

| Comparison w/ base models | Benchmarks
---|---|---
Dataset | GCN | R-GCN | GAT | R-GAT | SAGE |  R-SAGE | LAP | P-GCN | P-GAT | L-GCN | L-GAT
Cora | $72.7\atop\pm 1.35$ | $\bm{73.2}\atop\pm 1.20$ | $70.6\atop\pm 1.89$ | $72.0\atop\pm 1.77$ | $72.8\atop\pm 0.77$ |  $73.1\atop\pm 1.13$ | $72.8\atop\pm 1.22$ | $73.4\atop\pm 0.75$ | $72.8\atop\pm 1.34$ | $\bm{73.6}\atop\pm 1.15$ | $72.3\atop\pm 0.80$
Citeseer | $65.5\atop\pm 1.33$ | $\bm{66.4}\atop\pm 1.27$ | $63.4\atop\pm 2.57$ | $63.8\atop\pm 2.18$ | $64.6\atop\pm 1.31$ |  $65.0\atop\pm 1.15$ | $\bm{66.1}\atop\pm 0.98$ | $65.7\atop\pm 1.92$ | $63.0\atop\pm 1.24$ | $65.9\atop\pm 2.55$ | $63.5\atop\pm 1.80$
Pubmed | $73.2\atop\pm 1.57$ | $\bm{73.8}\atop\pm 0.95$ | $72.8\atop\pm 1.01$ | $73.8\atop\pm 1.15$ | $73.3\atop\pm 0.90$ |  $73.5\atop\pm 0.81$ | $73.5\atop\pm 1.52$ | $\bm{73.5}\atop\pm 1.19$ | $73.4\atop\pm 0.85$ | $73.4\atop\pm 1.49$ | $72.9\atop\pm 1.42$
PPI | $70.2\atop\pm 2.06$ | $71.3\atop\pm 0.51$ | $97.6\atop\pm 0.54$ | $\bm{98.2}\atop\pm 0.11$ | $72.6\atop\pm 1.37$ |  $\underline{72.6}\atop\pm 1.37$ | $70.8\atop\pm 0.70$ | $69.7\atop\pm 2.04$ | $98.0\atop\pm 0.23$ | $69.8\atop\pm 1.42$ | $\bm{98.2}\atop\pm 0.12$

###  IV-B Analysis and ablation studies

We analyze our model with R-GCN, which has significant gain as compared with
its base model (cf. Table I). The graph
${\mathcal{G}}_{0}=({\mathcal{V}}_{0},{\mathcal{E}}_{0})$ is the main
connected component of Cora. Specifically, we want to study whether R-GCN
indeed generates distributional graph signals with desired properties. For
smoothness (cf. Section III-A), we have interpreted earlier $\phi({\bm{O}})$
of the output features of ${\bm{O}}$ as weights of marginal distributions. For
the column $\phi({\bm{O}})_{:,1}$, we take the subvector indexed by
${\mathcal{V}}_{0}$, then normalize and denote it by ${\bm{x}}$ (analysis of
other columns are in Appendix D). We compute its Fourier transform
$\hat{{\bm{x}}}$ and inspect its high-frequency components. We show
$\hat{{\bm{x}}}$ for different epochs in Fig. 4. We also show the spectral
plots for the last epoch (epoch 200) of GCN and of the signal of (normalized)
ground truth labels. We see a clear shrinkage of high-frequency components for
R-GCN (epoch 200).

Figure 4: Spectral plots of (normalized) signals of probability weights and
ground truth labels.

For non-uniformity (cf. Section III-B), we collect in the set
${\mathbb{K}}_{\text{R-GCN}}$ (resp. ${\mathbb{K}}_{\text{GCN}}$) the
probability weights for all the label classes and nodes, i.e., $18956$ entries
of $\phi({\bm{O}})$ for epoch 200 of R-GCN (resp. GCN). Non-uniformity
suggests that ${\mathbb{K}}_{\text{R-GCN}}$ contains less values near the
average $1/7$ and more values near $1$, as compared with
${\mathbb{K}}_{\text{GCN}}$. To verify, we compare
$|{\mathbb{K}}_{\text{R-GCN}}\cap[1/7-\epsilon_{1},1/7+\epsilon_{1}]|$ with
$|{\mathbb{K}}_{\text{GCN}}\cap[1/7-\epsilon_{1},1/7+\epsilon_{1}]|$, and
$|{\mathbb{K}}_{\text{R-GCN}}\cap[1-\epsilon_{2},1]|$ with
$|{\mathbb{K}}_{\text{GCN}}\cap[1-\epsilon_{2},1]|$. The plots for different
choices of (small) $\epsilon_{1},\epsilon_{2}$ are shown in Fig. 5. The
results agree with our speculation.

Figure 5: Number of instances of output probability weights within a given
range.

Next, we conduct experiments on different R-GCN variants to validate the
effectiveness of the components of our model. Recall that the newly introduced
loss $\mathcal{L}_{0}$ has two components $\mathcal{L}_{1}$ and
$\mathcal{L}_{2}$, accounting for smoothness and non-uniformity of
distributional graph signals. We introduce modifications of R-GCN as follows:
(a) R1-GCN: we only retain $\mathcal{L}_{1}$ in the loss; (b) R2-GCN: we only
retain $\mathcal{L}_{2}$ in the loss, and (c) R3-GCN: we input ${\bm{O}}$
directly in the loss without applying $\phi$. In this case, the entries of
${\bm{O}}$ are not bounded, we only keep $\mathcal{L}_{1}$, which is positive
semi-definite, in the loss (cf. Section III-C). The results are shown in Table
IV. We see that R-GCN remains the most effective as compared with its
variants. This suggests that each component of R-GCN plays a useful role.

TABLE IV: Ablation study

| Cora | Citeseer | Pubmed
---|---|---|---
R-GCN | $\bm{83.4}\pm 0.24$ | $\bm{73.6}\pm 0.70$ | $\bm{79.9}\pm 0.40$
R1-GCN | $80.8\pm 0.69$ | $70.7\pm 0.67$ | $79.2\pm 0.56$
R2-GCN | $82.0\pm 0.65$ | $72.8\pm 0.92$ | $79.5\pm 0.61$
R3-GCN | $69.6\pm 1.75$ | $64.5\pm 1.79$ | $75.7\pm 0.57$

##  V Conclusion

In this paper, we introduce the notion of distributional graph signals and
total variations that measure the smoothness of such signals. Based on this
and the concept of non-uniformity, we propose a regularization scheme that can
be applied directly to enhance the performance of many GNN models. The method
is analogous to the regularization method of a smooth signal prior in GSP.

##  Appendix A Data statistics

In Table V, we provide statistics of datasets used in Section IV.

TABLE V: Dataset statistics

Dataset | Nodes | Edges | Classes | Features
---|---|---|---|---
Cora | 2708 | 5429 | 7 | 1433
Citeseer | 3327 | 4732 | 6 | 3703
Pubmed | 19717 | 44338 | 3 | 500
Photo | 7487 | 119043 | 8 | 745
CS | 18333 | 81894 | 15 | 6805
Disease | 1044 | 1043 | 2 | 1000
Airport | 3188 | 18631 | 4 | 4
PPI | Ave. 2373 | Ave. 34171 | 121 | 50

##  Appendix B Link prediction and graph classification

We follow Section III-D to apply the proposed regularization method to link
prediction and graph classification.

For link prediction [42, 33], we want to predict whether two nodes in a
network are likely to have a link. Suppose a base GNN model for link
prediction is given. We insert an R-module (Fig. 3 (b)) to the last node
feature matrix in the model pipeline. The loss of the R-module is added
directly to the original loss. We test on Cora and Airport datasets with base
models: GCN, GAT, HGCN, and GIL. The results are shown in Table VI. Except for
HGCN, we see that the regularization can enhance the performance of the base
models.

TABLE VI: Link prediction. The best performer is highlighted in blue.

| GCN | R-GCN | GAT | R-GAT | HGCN | R-HGCN | GIL | R-GIL
---|---|---|---|---|---|---|---|---
Cora | $88.1\atop\pm 1.08$ | $88.4\atop\pm 1.20$ | $88.8\atop\pm 1.41$ | $89.1\atop\pm 1.26$ | $93.3\atop\pm 0.43$ | $\underline{93.3}\atop\pm 0.43$ | $91.1\atop\pm 13.4$ | $\bm{94.6}\atop\pm 7.99$
Airport | $93.0\atop\pm 0.44$ | $93.2\atop\pm 0.37$ | $93.6\atop\pm 0.61$ | $\underline{93.6}\atop\pm 0.61$ | $97.6\atop\pm 0.13$ | $\underline{97.6}\atop\pm 0.13$ | $97.3\atop\pm 3.73$ | $\bm{97.7}\atop\pm 3.27$

We next consider graph classification. For such a task, we want to determine
the class label of each graph in a dataset containing multiple graphs. Graph
classification is closely related to the theoretical problem of graph
isomorphism test, and GIN [43] is a GNN model that explores such a connection.
We use GIN and the variant GIN2 with the learnable importance of the target
node compared to its neighbors, as base models. Similarly to link prediction,
we insert an R-module (Fig. 3 (b)) to the last node feature matrix in the
model pipeline. We test with bioinformatics datasets MUTAG and PTC [44],
following the protocol described in [43] and report the 10-fold cross
validation accuracy. Comparison results are shown in Table VII. The
regularization indeed works for both models.

TABLE VII: Graph classification. The best performer is highlighted in blue.

| GIN | R-GIN | GIN2 | R-GIN2
---|---|---|---|---
MUTAG | $87.7\pm 8.80$ | $\bm{89.9}\pm 6.40$ | $86.7\pm 6.77$ | $89.3\pm 7.98$
PTC | $64.6\pm 8.53$ | $64.8\pm 5.91$ | $64.8\pm 9.37$ | $\bm{66.3}\pm 9.63$

##  Appendix C Proofs of theoretical results

In this appendix, we discuss and prove a general result that implies Theorem
1. In addition, we also prove Lemma 1. We start with a computation of the
Wasserstein distance.

Suppose ${\mathbb{S}}=\\{s_{1},\ldots,s_{m}\\}$ is a finite discrete set and
$d$ is the discrete metric on ${\mathbb{S}}$. For
$\mu,\nu\in\mathcal{P}({\mathbb{S}})$, let $(\mu(s_{i}))_{1\leq i\leq n}$ and
$(\nu(s_{i}))_{1\leq i\leq n}$ be their respective probability weights.

######  Lemma 2.

$\displaystyle W(\mu,\nu)^{2}=\frac{1}{2}\sum_{1\leq i\leq
m}|\mu(s_{i})-\nu(s_{i})|.$

###### Proof.

The result is known [18]. For completeness, we provide a self-contained
elementary proof. Let $\gamma=\big{(}\gamma(s_{i},s_{j})\big{)}_{1\leq i,j\leq
m}$ be in $\Gamma(\mu,\nu)$. We have

$\displaystyle\begin{split}&\sum_{1\leq i\leq m}\sum_{1\leq j\leq
m}\gamma(s_{i},s_{j})d(s_{i},s_{j})^{2}\\\ &=\sum_{1\leq i\leq m}\sum_{1\leq
j\neq i\leq m}\gamma(s_{i},s_{j})\\\ &=\sum_{1\leq i\leq m}\left(\sum_{1\leq
j\leq m}\gamma(s_{i},s_{j})-\gamma(s_{i},s_{i})\right)\\\ &=\sum_{1\leq i\leq
m}\left(\mu(s_{i})-\gamma(s_{i},s_{i})\right)\\\ &\geq\sum_{1\leq i\leq
m}\left(\mu(s_{i})-\min(\mu(s_{i}),\nu(s_{i})\right).\end{split}$ (7)

As $W(\mu,\nu)^{2}$ is defined by taking the infimum of the left-hand side
over all $\gamma\in\Gamma(\mu,\nu)$, we have $W(\mu,\nu)^{2}\geq\sum_{1\leq
i\leq m}\big{(}\mu(s_{i})-\min(\mu(s_{i}),\nu(s_{i})\big{)}$. By the same
argument, we also have $W(\mu,\nu)^{2}\geq\sum_{1\leq i\leq
m}\big{(}\nu(s_{i})-\min(\mu(s_{i}),\nu(s_{i})\big{)}$. Summing up these two
inequalities, we have

$\displaystyle 2W(\mu,\nu)^{2}$ $\displaystyle\geq\sum_{1\leq i\leq
m}\big{(}\mu(s_{i})+\nu(s_{i})-2\min(\mu(s_{i}),\nu(s_{i})\big{)}$
$\displaystyle=\sum_{1\leq i\leq m}|\mu(s_{i})-\nu(s_{i})|.$

Therefore, to prove the lemma, it suffices to show that there is a $\gamma$
such that $\gamma(s_{i},s_{i})=\min\big{(}\mu(s_{i}),\nu(s_{i})\big{)}$. For
this, we prove a slightly more general claim: if non-negative numbers
$(x_{i})_{1\leq i\leq m}$ and $(y_{i})_{1\leq i\leq m}$ satisfy $\sum_{1\leq
i\leq m}x_{i}=\sum_{1\leq j\leq m}y_{i}=a$, then there are non-negative
$(z_{i,j})_{1\leq i,j\leq m}$ such that $\sum_{1\leq j\leq
m}z_{i,j}=x_{i},1\leq i\leq m$, $\sum_{1\leq i\leq m}z_{i,j}=y_{j},1\leq j\leq
m$, and $z_{i,i}=\min(x_{i},y_{i}),1\leq i\leq m$.

We prove this by induction on $m$. The case for $m=1$ is trivially true by
taking $z_{1,1}=x_{1}=y_{1}$. For $m\geq 2$, without loss of generality, we
assume that $x_{1}\geq y_{1}$ and $x_{2}\leq y_{2}$. Then we choose
$z_{1,1}=y_{1}$, $z_{2,2}=x_{2}$, $z_{1,j}=0,1<j\leq m$ and $z_{i,2}=0,1\leq
i\neq 2\leq m$. As a result, we form another two sequences of non-negative
numbers $x_{1}-y_{1},x_{3},\ldots,x_{m}$ and $y_{2}-x_{2},y_{3},\ldots,y_{m}$
with both summing to $a-x_{2}-y_{1}$. By the induction hypothesis, we are able
to find non-negative $(z^{\prime}_{i,j})_{1\leq i,j\leq m-1}$ for the two new
sequences of length $m-1$ each. It suffices to let
$z_{i,j}=z^{\prime}_{i-1,j-1}$ for $i>1$ or $j>2$ and
$z_{i,1}=z^{\prime}_{i-1,1}$ for $i>1$ (illustrated in Fig. 6). This proves
the claim and hence the lemma.

Figure 6: The relations between $z_{i,j}$ and $z^{\prime}_{i,j}$.

∎

To state and prove a general form of Theorem 1, we need to introduce a few
more notions. We fix marginal distributions $\mathscr{N}=\\{\mu_{i}:1\leq
i\leq n\\}$ with $\mu_{i}\in\mathcal{P}({\mathbb{S}})$. For any pair of nodes
$v_{i}$ and $v_{j}$ and $s\in{\mathbb{S}}$, define

$\displaystyle\rho_{i,j}(s)=\mu_{j}(s)/\mu_{i}(s)\text{ if
}\mu_{j}(s)\leq\mu_{i}(s)\text{ and }1\text{ otherwise.}$ (8)

More generally, if $P=(v_{i_{0}},\ldots,v_{i_{l}})$ is a directed path on
${\mathcal{G}}$ from $v_{i_{0}}$ to $v_{i_{l}}$, then

$\displaystyle\rho_{P}(s)=\prod_{0\leq j<l}\rho_{i_{j},i_{j+1}}(s).$

It is important to point out that $\rho_{P}$ can be computed directly as long
as $\mathscr{N}$ is given.

In the graph ${\mathcal{G}}$, suppose ${\mathcal{H}}$ is a spanning tree and
$v_{0}$ is a fixed (root) node. Let ${\mathcal{E}}_{{\mathcal{H}}}$ be the
edge set of ${\mathcal{H}}$ and
${\mathcal{E}}^{\prime}={\mathcal{E}}\backslash{\mathcal{E}}_{{\mathcal{H}}}$.
For each edge $e=(v_{i},v_{j})$, let $P_{i}$ (resp. $P_{j}$) be the unique
path on ${\mathcal{H}}$ connecting $v_{0}$ and $v_{i}$ (resp. $v_{j}$).
Moreover, $v_{0}$ is an endpoint of $P_{i}\cap P_{j}$, and let $v_{k}$ be the
other endpoint of $P_{i}\cap P_{j}$. Denote by $Q_{i}$ (resp. $Q_{j}$) be the
direct path (on ${\mathcal{H}}$) from $v_{k}$ to $v_{i}$ (resp. $v_{j}$) (see
Fig. 7). We introduce

$\displaystyle\mathfrak{t}_{i,j}(s)=\mu_{i}(s)+\mu_{j}(s)-2\mu_{k}(s)\rho_{Q_{i}}(s)\rho_{Q_{j}}(s).$
(9)

######  Definition 6.

Define

$\displaystyle\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})=\sum_{s\in{\mathbb{S}}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}\mathfrak{t}_{i,j}(s).$

Figure 7: An example of paths $Q_{i}$ and $Q_{j}$.

We can compute the special case where ${\mathcal{G}}=\mathcal{H}$ is a tree.
Notice that for an edge $e=(v_{i},v_{j})\in{\mathcal{E}}_{\mathcal{H}}$
directed from $v_{i}$ to $v_{j}$, then $v_{k}=v_{i}$ and $Q_{i}=\\{v_{i}\\}$,
$Q_{j}=e$ and $\rho_{Q_{i}}(s)=1$. If $\mu_{i}(s)\geq\mu_{j}(s)$, then
$2\mu_{k}(s)\rho_{Q_{i}}(s)\rho_{Q_{j}}(s)=2\mu_{i}(s)\cdot\mu_{j}(s)/\mu_{i}(s)=2\mu_{j}(s)$.
Hence, we have

$\displaystyle\mu_{i}(s)$
$\displaystyle+\mu_{j}(s)-2\mu_{k}(s)\rho_{Q_{i}}(s)\rho_{Q_{j}}(s)$
$\displaystyle=\mu_{i}(s)-\mu_{j}(s)=|\mu_{i}(s)-\mu_{j}(s)|.$

The case $\mu_{i}(s)<\mu_{j}(s)$ is similar, and in summary

$\begin{split}\mathfrak{t}_{i,j}(s)&=\mu_{i}(s)+\mu_{j}(s)-2\mu_{k}(s)\rho_{Q_{i}}(s)\rho_{Q_{j}}(s)\\\
&=|\mu_{i}(s)-\mu_{j}(s)|\end{split}$ (10)

Therefore, for any $v_{0}$, we have

$\displaystyle\begin{split}\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})&=\sum_{s\in{\mathbb{S}}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}|\mu_{i}(s)-\mu_{j}(s)|\\\
&=\sum_{s\in{\mathbb{S}}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}|\bm{\mu}_{s}(i)-\bm{\mu}_{s}(j)|=\mathcal{T}_{G,1}(\mathscr{N}),\end{split}$
(11)

where $\bm{\mu}_{s}(i)$ is same as $\mu_{i}(s)$ (cf. Section III-A).

Following the notations in Section III-A, we have the following generalization
of Theorem 1.

######  Theorem 2.

For any spanning tree ${\mathcal{H}}$ of ${\mathcal{G}}$ and root node
$v_{0}$, we have

$\displaystyle\mathcal{T}_{G,2}(\mathscr{N})\leq\mathcal{T}_{G,1}(\mathscr{N})\leq
2\mathcal{T}_{G}(\mathscr{N})\leq\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N}).$

Before proving the result, we remark that although
$\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})$ is defined in a convoluted
way, it can however be computed directly given $\mathscr{N},\mathcal{H}$ and
$v_{0}$. Therefore, the result gives computable upper and lower bounds of
$\mathcal{T}_{G}(\mathscr{N})$.

###### Proof.

As $|\bm{\mu}_{s}(i)-\bm{\mu}_{s}(j)|\leq 1$, it is trivially true that
$\mathcal{T}_{G,2}(\mathscr{N})\leq\mathcal{T}_{G,1}(\mathscr{N})$.

To show $\mathcal{T}_{G,1}(\mathscr{N})\leq 2\mathcal{T}_{G}(\mathscr{N})$, we
first claim that
$\mathcal{T}_{G}(\mathscr{N})=\inf_{\mu\in\Gamma(\mathscr{N})}\mathcal{T}_{G}(\mu)$
is achieved for some $\mu_{0}\in\Gamma(\mathscr{N})$. The map
$\alpha:\Gamma(\mathscr{N})\to\mathbb{R},\mu\mapsto\mathcal{T}_{G}(\mu)$ is
continuous. On the other hand, $\Gamma(\mathscr{N})$ is a compact subset of a
Euclidean space. This is because $\mathcal{P}({\mathbb{S}}^{n})$ is a bounded
subset of $\mathbb{R}^{m^{n}}$ with the components corresponding to weights of
the joint distribution. Moreover, $\Gamma(\mathscr{N})$ is closed because the
condition to have the prescribed marginal distributions is a set of linear
conditions. By the extreme value theorem, $\inf\alpha$ is achieved for some
$\mu_{0}\in\Gamma(\mathscr{N})$. For each edge
$(v_{i},v_{j})\in{\mathcal{E}}$, let $\mu_{0,i,j}$ be the marginal
distribution of the pair $(v_{i},v_{j})$. The marginals of $\mu_{0,i,j}$ are
$\mu_{i}$ and $\mu_{j}$ at $v_{i}$ and $v_{j}$ respectively. We have

$\displaystyle\mathcal{T}_{G}(\mathscr{N})$
$\displaystyle=\mathcal{T}_{G}(\mu_{0})=\int\sum_{(v_{i},v_{j})\in{\mathcal{E}}}d({\bm{x}}(i),{\bm{x}}(j))^{2}\,\mathrm{d}\mu_{0}({\bm{x}})$
$\displaystyle=\sum_{(v_{i},v_{j})\in{\mathcal{E}}}\int
d({\bm{x}}(i),{\bm{x}}(j))^{2}\,\mathrm{d}\mu_{0}({\bm{x}})$
$\displaystyle=\sum_{(v_{i},v_{j})\in{\mathcal{E}}}\int
d({\bm{y}}(i),{\bm{y}}(j))^{2}\,\mathrm{d}\mu_{0,i,j}({\bm{y}})$
$\displaystyle\stackrel{{\scriptstyle\text{Def.}~{}\ref{defn:fmm}}}{{\geq}}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}W(\mu_{i},\mu_{j})^{2}$
$\displaystyle\stackrel{{\scriptstyle\text{Lem.}~{}\ref{lem:wmn}}}{{=}}\frac{1}{2}\sum_{(v_{i},v_{j})\in{\mathcal{E}}}\sum_{s\in{\mathbb{S}}}|\bm{\mu}_{s}(i)-\bm{\mu}_{s}(j)|=\frac{1}{2}\mathcal{T}_{G,1}(\mathscr{N}).$

We now prove
$2\mathcal{T}_{G}(\mathscr{N})\leq\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})$.
Given a spanning tree $\mathcal{H}$ and node $v_{0}$, we construct
$\mu_{\mathcal{H},v_{0}}\in\Gamma(\mathscr{N})$ using ideas from the theory of
Bayesian networks [45] as follows. We make $\mathcal{H}$ directed by requiring
that each edge is pointed away from the (root) node $v_{0}$. As a consequence,
each node has at most $1$ incoming edge. Let ${\bm{x}}=({\bm{x}}_{i})_{1\leq
i\leq n}$ be the random vector with ${\bm{x}}_{i}$ the (random) label at
$v_{i}$. For each directed edge $(v_{i},v_{j})$ in $\mathcal{H}$, let
$\mu_{i,j}\in\Gamma(\mu_{i},\mu_{j})$ be a distribution that realizes
$W(\mu_{i},\mu_{j})$, which give a conditional probability weights

$\displaystyle{\bm{p}}_{i,j}$
$\displaystyle=\\{p_{i,j}(s,s^{\prime}):s,s^{\prime}\in{\mathbb{S}}\\}$
$\displaystyle=\\{p({\bm{x}}_{j}=s^{\prime}\mid{\bm{x}}_{i}=s):s,s^{\prime}\in{\mathbb{S}}\\}.$

By [45] Section 8.1 (8.5), there is a distribution
$\mu_{\mathcal{H},v_{0}}\in\Gamma(\mathscr{N})$ such that its marginal for
each edge $(v_{i},v_{j})\in{\mathcal{E}}_{\mathcal{H}}$ is $\mu_{i,j}$. For
each edge $(v_{i},v_{j})\in{\mathcal{E}}^{\prime}$, let $\mu_{i,j}^{\prime}$
be the marginal of $\mu_{\mathcal{H},v_{0}}$ to the pair $(v_{i},v_{j})$. As
$\mu_{i,j}$ realizes $W(\mu_{i},\mu_{j})$, by Lemma 2,
$\mu_{i,j}(s,s)=\min(\mu_{i}(s),\mu_{j}(s))$ for each $s\in{\mathbb{S}}$. In
particular, $p_{i,j}(s,s)=\rho_{i,j}(s)$ (cf. (8)). More generally, if $P$ is
a directed path in $\mathcal{H}$ from $v_{i}$ to $v_{j}$, then the following
inequality holds

$\displaystyle p_{i,j}(s,s)\geq\rho_{P}(s).$ (12)

According to the definition, we have
$\mathcal{T}_{G}(\mathscr{N})\leq\mathcal{T}_{G}(\mu_{\mathcal{H},v_{0}})=\mathcal{S}_{{\mathcal{E}}_{\mathcal{H}}}+\mathcal{S}_{{\mathcal{E}}^{\prime}}$.
The summand

$\displaystyle\mathcal{S}_{{\mathcal{E}}_{\mathcal{H}}}=\sum_{(v_{i},v_{j})\in{\mathcal{E}}_{\mathcal{H}}}\int
d({\bm{y}}(i),{\bm{y}}(j))^{2}\,\mathrm{d}\mu_{i,j}({\bm{y}})$

is the summation over the edges of ${\mathcal{E}}_{\mathcal{H}}$, while

$\displaystyle\mathcal{S}_{{\mathcal{E}}^{\prime}}=\sum_{(v_{i},v_{j})\in{\mathcal{E}}^{\prime}}\int
d({\bm{y}}(i),{\bm{y}}(j))^{2}\,\mathrm{d}\mu_{i,j}^{\prime}({\bm{y}})$

is the summand over ${\mathcal{E}}^{\prime}$. For
$\mathcal{S}_{{\mathcal{E}}_{\mathcal{H}}}$, we have seen that for each edge
$(v_{i},v_{j})\in{\mathcal{E}}_{\mathcal{H}}$

$\displaystyle\int
d({\bm{y}}(i),{\bm{y}}(j))^{2}\,\mathrm{d}\mu_{i,j}({\bm{y}})$
$\displaystyle\stackrel{{\scriptstyle\text{Lem.}~{}\ref{lem:wmn}}}{{=}}$
$\displaystyle\frac{1}{2}\sum_{s\in{\mathbb{S}}}|\mu_{i}(s)-\mu_{j}(s)|$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:mmm})}}{{=}}$
$\displaystyle\frac{1}{2}\sum_{s\in{\mathbb{S}}}\big{(}\mu_{i}(s)+\mu_{j}(s)-2\mu_{k}(s)\rho_{Q_{i}}(s)\rho_{Q_{j}}(s)\big{)},$

where the right-hand-side is the term $\frac{1}{2}\mathfrak{t}_{i,j}(s)$ (cf.
(9)) that corresponds to $(v_{i},v_{j})$ in
$\frac{1}{2}\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})$.

Consider $(v_{i},v_{j})\in{\mathcal{E}}^{\prime}$ and we first notice the
following identity:

$\displaystyle 2$ $\displaystyle\int
d({\bm{y}}(i),{\bm{y}}(j))^{2}\,\mathrm{d}\mu^{\prime}_{i,j}({\bm{y}})$
$\displaystyle\stackrel{{\scriptstyle(\ref{eq:sum})}}{{=}}\sum_{s\in{\mathbb{S}}}\mu_{i}(s)+\mu_{j}(s)-2\mu^{\prime}_{i,j}(s,s).$

To show the summand is bounded by $\mathfrak{t}_{i,j}(s)$, it suffices to show
$\mu_{i,j}^{\prime}(s,s)$ is bounded below by
$\rho_{Q_{i}}(S)\rho_{Q_{j}}(s)\mu_{k}(s)$, with $v_{k}$ and paths
$Q_{i},Q_{j}$ be as in (9). We estimate using the construction of
$\mu_{i,j}^{\prime}$ based on the method of Bayesian network (on directed
$\mathcal{H}$) as follows:

$\displaystyle\mu^{\prime}_{i,j}(s,s)$ $\displaystyle\geq
p({\bm{x}}_{i}=s,{\bm{x}}_{j}=s,{\bm{x}}_{k}=s)$
$\displaystyle=p({\bm{x}}_{i}=s,{\bm{x}}_{j}=s\mid{\bm{x}}_{k}=s)\mu_{k}(s)$
$\displaystyle=p_{k,i}(s,s)p_{k,j}(s,s)\mu_{k}(s)$
$\displaystyle\geq\rho_{Q_{i}}(s)\rho_{Q_{j}}(s)\mu_{k}(s).$

For the last equality, we use the fact that $Q_{i}\cap Q_{j}=\\{v_{k}\\}$; and
hence ${\bm{x}}_{i},{\bm{x}}_{j}$ are independent given ${\bm{x}}_{k}$ by the
Bayesian network construction. The last line follows from (12). Consequently,
the inequality that
$2\mathcal{T}_{G}(\mathscr{N})\leq\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})$
follows. ∎

If we examine the formula of
$\mathcal{T}_{G,{\mathcal{H}},v_{0}}(\mathscr{N})$ (by changing summation
order), it can be decomposed into two parts
$\sum_{(v_{i},v_{j})\in{\mathcal{E}}_{\mathcal{H}}}$ and
$\sum_{(v_{i},v_{j})\in{\mathcal{E}}^{\prime}}$. The former is essentially
$\mathcal{T}_{G,1}$ of $\mathscr{N}$ on the tree $\mathcal{H}$. This is the
reason that Theorem 2 implies Theorem 1. Formally we have the following.

###### Proof of Theorem 1..

For $G=(V,E)$ of size $n$, recall $K_{n}$ is the complete graph on $V$ and the
complement $\overline{G}$ of $G$ is obtained from $K_{n}$ by removing edges in
$E$. Let $\omega(G)$ be the clique number of $\overline{G}$, i.e., the number
of vertices in a largest clique of $\overline{G}$. Let $c_{1}=n-\omega(G)$,
$c_{3}=\sqrt{|{\mathbb{S}}|n}$ and $c_{2}=c_{1}c_{3}$. Notice that
$c_{1},c_{2},c_{3}$ are independent of $\mathscr{N}$.

It suffices to prove the following inequalities

  * •

$\mathcal{T}_{G,1}(\mathscr{N})\leq 2\mathcal{T}^{c}_{G}(\mathscr{N})$;

  * •

$2\mathcal{T}^{c}_{G}(\mathscr{N})\leq c_{1}\mathcal{T}_{G,1}(\mathscr{N})$;

  * •

$\mathcal{T}_{G,1}(\mathscr{N})\leq
c_{3}\mathcal{T}_{G,2}(\mathscr{N})^{\frac{1}{2}}$.

Let ${\mathbb{T}}=\\{T_{1},\ldots,T_{k}\\}$ be any cover of $G$ by spanning
trees. We have $\mathcal{T}_{G,1}(\mathscr{N})\leq\sum_{1\leq i\leq
k}\mathcal{T}_{G,1}(T_{i})=\sum_{1\leq i\leq
k}2\mathcal{T}_{T_{i}}(\mathscr{N})$ by Theorem 2. As the inequality is true
for any ${\mathbb{T}}$, we have $\mathcal{T}_{G,1}(\mathscr{N})\leq
2\mathcal{T}^{c}_{G}(\mathscr{N})$.

For the second inequality, by the theory of tree cover numbers of graphs [46,
Prop. 2], there is a ${\mathbb{T}}=\\{T_{1},\ldots,T_{c_{1}}\\}$ that covers
$G$. Therefore, by Theorem 2,
$2\mathcal{T}^{c}_{G}(\mathscr{N})\leq\sum_{1\leq i\leq
c_{1}}\mathcal{T}_{T_{i}}(\mathscr{N})\leq
c_{1}\mathcal{T}_{G,1}(\mathscr{N})$.

The last inequality follows from the Cauchy–Schwarz inequality. ∎

We end this section by proving Lemma 1.

###### Proof of Lemma 1..

Using Lemma 2, we estimate

$\displaystyle 2\sum_{1\leq i\leq
n}a_{i}W(\mu_{{\bm{X}},i},\mathcal{U}({\mathbb{S}}))^{2}$ $\displaystyle=$
$\displaystyle\sum_{1\leq i\leq n}a_{i}\sum_{1\leq j\leq
m}\left\lvert{\bm{X}}_{i,j}-\frac{1}{|{\mathbb{S}}|}\right\rvert$
$\displaystyle\leq$ $\displaystyle\sum_{1\leq i\leq n}a_{i}\sum_{1\leq j\leq
m}({\bm{X}}_{i,j}-\frac{1}{|{\mathbb{S}}|})^{2}$ $\displaystyle=$
$\displaystyle\sum_{1\leq i\leq n}a_{i}\sum_{1\leq j\leq
m}{\bm{X}}_{i,j}^{2}-2\sum_{1\leq i\leq n}a_{i}\sum_{1\leq j\leq
m}\frac{{\bm{X}}_{i,j}}{|{\mathbb{S}}|}$ $\displaystyle+\sum_{1\leq i\leq
n}a_{i}\sum_{1\leq j\leq m}\frac{1}{|{\mathbb{S}}|^{2}}$ $\displaystyle=$
$\displaystyle\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})-2\sum_{1\leq
i\leq n}a_{i}\frac{1}{|{\mathbb{S}}|}+\sum_{1\leq i\leq
n}a_{i}\frac{1}{|{\mathbb{S}}|}$ $\displaystyle=$
$\displaystyle\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})-\frac{1}{|{\mathbb{S}}|}\operatorname{Tr}({\bm{D}}).$

Therefore, $2\sum_{1\leq i\leq
n}a_{i}W(\mu_{{\bm{X}},i},\mathcal{U}({\mathbb{S}}))^{2}\leq\operatorname{Tr}({\bm{X}}^{\top}{\bm{D}}{\bm{X}})+C$
with $C=-\frac{1}{|{\mathbb{S}}|}\operatorname{Tr}({\bm{D}})$, which is
independent of ${\bm{X}}$.

The last statement of the lemma follows from the simple fact that for each
$1\leq i\leq n$, we have

$\displaystyle\sum_{1\leq j\leq m}\frac{1}{|{\mathbb{S}}|^{2}}\leq\sum_{1\leq
j\leq m}{\bm{X}}_{i,j}^{2}\leq 1.$

∎

##  Appendix D Additional plots for model analysis

We supplement Fig. 4 by showing spectral plots of signals of probability
weights $\phi({\bm{O}})_{:,i},2\leq i\leq 7$ for the Cora dataset. The index
$i$ corresponds to the $i$-th label class. From the plots, we observe that
during the training of R-GCN, the high-frequency components indeed shrink for
all the label classes. Compared with GCN, the last epoch of R-GCN has smaller
high-frequency components for the $3$rd and $5$th label classes

Figure 8: Spectral plots of (normalized) signals of probability weights for
the $2$nd label class.

Figure 9: Spectral plots of (normalized) signals of probability weights for
the $3$rd label class.

Figure 10: Spectral plots of (normalized) signals of probability weights for
the $4$th label class.

Figure 11: Spectral plots of (normalized) signals of probability weights for
the $5$th label class.

Figure 12: Spectral plots of (normalized) signals of probability weights for
the $6$th label class.

Figure 13: Spectral plots of (normalized) signals of probability weights for
the $7$th label class.

## References

  * [1] T. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _International Conference on Learning Representations_ , 2017.
  * [2] M. Defferrard, X. Bresson, and P. Vandergheynst, “Convolutional neural networks on graphs with fast localized spectral filtering,” in _Advances in Neural Information Processing Systems_ , 2016, p. 3844–3852.
  * [3] P. Velic̄ković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph attention networks,” in _International Conference on Learning Representations_ , 2018.
  * [4] I. Chami, Z. Ying, C. Ré, and J. Leskovec, “Hyperbolic graph convolutional neural networks,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 4869–4880.
  * [5] W. L. Hamilton, R. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in _Advances in Neural Information Processing Systems_ , 2017, p. 1025–1035.
  * [6] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst, “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,” _IEEE Signal Process. Mag._ , vol. 30, no. 3, pp. 83–98, 2013.
  * [7] A. Ortega, P. Frossard, J. Kovačević, J. M. F. Moura, and P. Vandergheynst, “Graph signal processing: Overview, challenges, and applications,” _Proc. IEEE_ , vol. 106, no. 5, pp. 808–828, 2018.
  * [8] F. Ji and W. P. Tay, “A Hilbert space theory of generalized graph signal processing,” _IEEE Trans. Signal Process._ , vol. 67, no. 24, pp. 6188 – 6203, 2019.
  * [9] S. Narang, A. Gadde, and A. Ortega, “Signal processing techniques for interpolation in graph structured data,” in _International Conference on Acoustics, Speech, and Signal Processing_ , 2013.
  * [10] S. Chen, A. Sandryhaila, J. M. F. Moura, and J. Kovac̄ević, “Signal recovery on graphs: Variation minimization,” _IEEE Trans. Signal Process._ , vol. 63, no. 17, pp. 4609–4624, 2015.
  * [11] V. Kalofolias, “How to learn a graph from smooth signals,” in _International Conference on Artificial Intelligence and Statistics_ , 2016, pp. 920–929.
  * [12] H. Yang, K. Ma, and J. Cheng, “Rethinking graph regularization for graph neural networks,” in _AAAI Conference on Artificial Intelligence_ , 2021, pp. 4573–4581.
  * [13] F. Feng, X. He, J. Tang, and T. Chua, “Graph adversarial training: Dynamically regularizing based on graph structure,” _IEEE Trans. Knowl. Data Eng._ , vol. 33, no. 6, pp. 2493–2504, 2021.
  * [14] Y. Deng, Z. Dong and J. Zhu, “Batch virtual adversarial training for graph convolutional networks,” _arXiv:1902.09192_ , 2019.
  * [15] O. Stretcu, K. Viswanathan, D. Movshovitz-Attias, E. Platanios, S. Ravi, and A. Tomkins, “Graph agreement models for semi-supervised learning,” in _Advances in Neural Information Processing Systems_ , 2019, p. 8713–8723.
  * [16] P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Galligher, and T. Eliassi-Rad, “Collective classification in network data,” _AI Mag._ , vol. 29, no. 3, p. 93, 2008.
  * [17] F. Ji, X. Jian, and W. P. Tay, “On distributional graph signals,” _arXiv preprint arXiv:2302.11104_ , 2023.
  * [18] C. Villani, _Optimal Transport, Old and New_. Springer, 2009.
  * [19] D. Shah and T. Zaman, “Rumors in a network: Who’s the culprit?” _IEEE Trans. Inf. Theory_ , vol. 57, no. 8, pp. 5163–5181, 2011.
  * [20] F. Ji, W. Tang, and W. P. Tay, “On the properties of gromov matrices and their applications in network inference,” _IEEE Trans. Signal Process._ , vol. 67, no. 10, pp. 2624–2638, 2019.
  * [21] X. Ma, H. Chen, and G. Song, “LEReg: Empower graph neural networks with local energy regularization,” in _Conference on Information and Knowledge Management_ , 2021.
  * [22] X. Zhu, Z. Ghahramani, and J. Lafferty, “Semi-supervised learning using Gaussian fields and harmonic functions,” in _International Conference on Machine Learning_ , 2003, p. 912–919.
  * [23] D. Zhou, O. Bousquet, T. Lal, J. Weston, and B. Schölkopf, “Learning with local and global consistency,” in _Advances in Neural Information Processing Systems_ , 2004, p. 321–328.
  * [24] R. Ando and T. Zhang, “Learning on graph with Laplacian regularization,” in _Advances in Neural Information Processing Systems_ , 2007, pp. 25–32.
  * [25] H. NT and T. Maehara, “Revisiting graph neural networks: All we have is low-pass filters,” _arXiv:1905.09550_ , 2019.
  * [26] K. Oono and T. Suzuki, “Graph neural networks exponentially lose expressive power for node classification,” in _International Conference on Learning Representations_ , 2020.
  * [27] L. Zhao and L. Akoglu, “PairNorm: Tackling oversmoothing in GNNs,” in _International Conference on Learning Representations_ , 2020.
  * [28] D. Chen, Y. Lin, W. Li, P. Li, J. Zhou, and X. Sun, “Measuring and relieving the over-smoothing problem for graph neural networks from the topological view,” in _AAAI Conference on Artificial Intelligence_ , 2020, pp. 3438–3445.
  * [29] Y. Rong, W. Huang, T. Xu, and J. Huang, “DropEdge: Towards deep graph convolutional networks on node classification,” in _International Conference on Learning Representations_ , 2020.
  * [30] G. Li, M. Müller, A. Thabet, and B. Ghanem, “DeepGCNs: Can GCNs go as deep as CNNs?” in _International Conference on Computer Vision_ , 2019\.
  * [31] S. Luan, M. Zhao, X. Chang, and D. Precup, “Break the ceiling: Stronger multi-scale deep graph convolutional networks,” in _Advances in Neural Information Processing Systems_ , 2019, pp. 10 943–10 953.
  * [32] G. Namata, B. London, L. Getoor, and B. Huang, “Query-driven active surveying for collective classification,” in _Mining and Learning with Graphs_ , 2012\.
  * [33] M. Zhang and Y. Chen, “Link prediction based on graph neural networks,” in _Advances in Neural Information Processing Systems_ , 2018, p. 5171–5181.
  * [34] O. Shchur, M. Mumme, A. Bojchevski, and S. Günnemann, “Pitfalls of graph neural network evaluation,” in _Advances in Neural Information Processing Systems_ , 2018.
  * [35] M. Fey and J. Lenssen, “Fast graph representation learning with PyTorch Geometric,” in _International Conference on Learning Representations_ , 2019.
  * [36] D. Szklarczyk, J. Morris, H. Cook, M. Kuhn, S. Wyder, M. Simonovic, A. Santos, N. Doncheva, A. Roth, P. Bork, and et al., “The string database in 2017: quality-controlled protein–protein association networks, made broadly accessible,” _Nucleic acids research_ , 2016.
  * [37] T. Rusch, B. Chamberlain, J. Rowbottom, S. Mishra, and M. Bronstein, “Graph-coupled oscillator networks,” in _International Conference on Machine Learning_ , 2022.
  * [38] M. Bridson and A. Haefliger, _Metric Spaces of Non-positive Curvature_. Springer, 1999.
  * [39] C. Gulcehre, M. Denil, M. Malinowski, A. Razavi, R. Pascanu, K. M. Hermann, P. Battaglia, V. Bapst, D. Raposo, and A. Santoro, “Hyperbolic attention networks,” in _International Conference on Learning Representations_ , 2019\.
  * [40] S. Zhu, S. Pan, C. Zhou, J. Wu, Y. Cao, and B. Wang, “Graph geometry interaction learning,” in _Advances in Neural Information Processing Systems_ , 2020, p. 7548–7558.
  * [41] P. Mishra, A. Piktus, G. Goossen, and F. Silvestri, “Node masking: Making graph neural networks generalize and scale better,” _arXiv:2001.07524v4_ , 2021.
  * [42] D. Liben-Nowell and J. Kleinberg, “The link-prediction problem for social networks,” _JASIST_ , vol. 58, no. 7, pp. 1019–1031, 2007.
  * [43] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful are graph neural networks?” in _International Conference on Learning Representations_ , 2019\.
  * [44] P. Yanardag and S. Vishwanathan, “Deep graph kernels,” in _SIGKDD Conference on Knowledge Discovery and Data Mining_ , 2015.
  * [45] C. Bishop, _Pattern Recognition and Machine Learning_. Springer, 2006.
  * [46] C. Bozeman, M. Catral, B. Cook, O. González, and C. Reinhart, “On the tree cover number of a graph,” _Involve: A Journal of Mathematics_ , vol. 10, no. 5, pp. 767 – 779, 2017.