# Understanding Non-linearity in Graph Neural Networks from the Perspective of
Bayesian Inference

Rongzhe Wei1, Haoteng Yin1, Junteng Jia2, Austin R. Benson2, Pan Li1
1 Department of Computer Science, Purdue University
2 Department of Computer Science, Cornell University
{wei397, yinht<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Graph neural networks (GNNs) have shown superiority in many prediction tasks
over graphs due to their impressive capability of capturing nonlinear
relations in graph-structured data. However, for node classification tasks,
often, only marginal improvement of GNNs over their linear counterparts has
been observed. Previous works provide very few understandings of this
phenomenon. In this work, we resort to Bayesian learning to deeply investigate
the functions of non-linearity in GNNs for node classification tasks. Given a
graph generated from the statistical model CSBM, we observe that the max-a-
posterior estimation of a node label given its own and neighbors’ attributes
consists of two types of non-linearity, a possibly non-linear transformation
of node attributes and a ReLU-activated feature aggregation from neighbors.
The latter surprisingly matches the type of non-linearity used in many GNN
models. By further imposing a Gaussian assumption on node attributes, we prove
that the superiority of those ReLU activations is only significant when the
node attributes are far more informative than the graph structure, which
nicely matches many previous empirical observations. A similar argument can be
achieved when there is a distribution shift of node attributes between the
training and testing datasets. Finally, we verify our theory on both synthetic
and real-world networks. Our code is available at https://github.com/Graph-
COM/Bayesian_inference_based_GNN.git.

##  1 Introduction

Learning on graphs (LoG) has been widely used in the applications with graph-
structured data [zhu2005semi, GRLbook]. Node classification, as one of the
most crucial tasks in LoG, asks to predict the labels of nodes in a graph,
which has been used in many applications such as community detection
[fortunato2010community, lancichinetti2009community, chen2020supervised,
liu2021deep], anomaly detection [ma2021comprehensive, wang2021bipartite],
biological pathway analysis [aittokallio2006graph, scott2006efficient] and so
on.

Recently, graph neural networks (GNNs) have become the de-facto standard used
in many LoG tasks due to their super empirical performance
[hamilton2017inductive, kipf2016semi]. Researchers often attribute such
success to non-linearity in GNNs which associates them with great expressive
power [xu2018powerful, morris2019weisfeiler]. GNNs can approximate a wide
range of functions defined over graphs [keriven2019universal,
chen2019equivalence, azizian2021expressive] and thus excel in predicting,
e.g., the free energies of molecules [gilmer2017neural], which are by nature
non-linear solutions of some quantum-mechanical equations. However, for node
classification tasks, many studies have shown that non-linearity to control
the exchange of features among neighbors seems not that crucial. For example,
many works use linear propagation of node attributes over graphs
[wu2019simplifying, he2020lightgcn], and others recommend adding non-linearity
while only to the transformation of initial node attributes
[klicpera2018predict, DBLP:conf/nips/KlicperaWG19, huang2020combining]. Both
cases achieve comparable or even better performance than other models with
complex nonlinear propagation, such as using neighbor-attention mechanism
[velivckovic2018graph]. Recently, even in the complicated heterophilic setting
where nodes with same labels are not directly connected, linear propagation
still achieves competitive performance [chien2021adaptive, wang2022powerful],
compared with the models with nonlinear and deep architectures [zhu2020beyond,
zhu2021graph].

Although empirical studies on GNNs are extensive till now and many practical
observations as above have been made, there have been very few works
attempting to characterize GNNs in theory, especially to understand the effect
of non-linearity by comparing with the linear counterparts for node
classification tasks. The only work on this topic to the best of our knowledge
still focuses on comparing the expressive power of the two methods to
distinguish nodes with different local structures [chen2020graph]. However,
the achieved statement that non-linear propagation improves expressiveness may
not necessarily reveal the above phenomenon that non-linear and linear methods
have close empirical performance while with subtle difference. Moreover, more
expressiveness is often at the cost of model generalization and thus may not
necessarily yield more accurate prediction [cong2021provable,
li2022expressive].

In this work, we expect to give a more precise characterization of the values
of non-linearity in GNNs from a statistical perspective, based on Bayesian
inference specifically. We resort to contextual stochastic block models (CSBM)
[binkiewicz2017covariate, deshpande2018contextual]. We make a significant
observation that given a graph generated by CSBM, the max-a-posterior (MAP)
estimation of a node label given its own and neighbors’ features surprisingly
corresponds to a graph convolution layer with _ReLU_ as the activation
combined with an initial node-attribute transformation. Such a transformation
of node attributes is generally nonlinear unless they are generated from the
natural exponential family [morris1982natural]. Since the MAP estimator is
known to be Bayesian optimal [theodoridis2015machine], the above observation
means that ReLU-based propagation has the potential to outperform linear
propagation. To precisely characterize such benefit, we further assume that
the node attributes are generated from a label-conditioned Gaussian model, and
analyze and compare the node mis-classification errors of linear and nonlinear
models. We have achieved the following conclusions (note that we only provide
informal statements here and the formal statements are left in the theorems).

  * •

When the node attributes are less informative compared to the structural
information, non-linear propagation and linear propagation have almost the
same mis-classification error (case I in Thm. 2).

  * •

When the node attributes are more informative, non-linear propagation shows
advantages. The mis-classification error of non-linear propagation can be
significantly smaller than that of linear propagation with sufficiently
informative node attributes (case II in Thm. 2).

  * •

When there is a distribution shift of the node attributes between the training
and testing datasets, non-linearity provides better transferability in the
regime of informative node attributes (Thm. 3).

Given that practical node attributes are often not that informative, the
advantages of non-linear propagation over linear propagation for node
classification is limited albeit observable. Our analysis and conclusion apply
to both homophilic and heterophilic settings, i.e., when nodes with same
labels tend to be connected (homophily) or disconnected (heterophily),
respectively [chien2021adaptive, zhu2020beyond, zhu2021graph,
suresh2021breaking, ma2021homophily].

Figure 1: Averaged one-vs-all Classification Accuracy on Citation Networks of
Nonlinear Models v.s. Linear Models. Node attributes in or out of the one
class are generated from Gaussian distributions $\mathcal{N}(\mu,\frac{I}{m})$
and $\mathcal{N}(\nu,\frac{I}{m})$, $\mu,\nu\in\mathbb{R}^{m}$, respectively.
The detailed settings are introduced in Sec. LABEL:sec:exp-real.

Extensive evaluation on both synthetic and real datasets demonstrates our
theory. Specifically, the node mis-classification errors of three citation
networks with different levels of attributed information (Gaussian attributes)
are shown in Fig. 1, which precisely matches the above conclusions.

###  1.1 More Related Works

GNNs have achieved great empirical success while theoretical understanding of
GNNs, their non-linearity in particular, is still limited. There are many
works studying the expressive power of GNNs [maron2019provably,
balcilar2021breaking, azizian2020expressive, murphy2019relational,
sato2021random, abboud2021surprising, bodnar2021weisfeiler,
chen2019equivalence, li2020distance, loukas2020graph, vignac2020building,
zhang2021nested], while they often assume arbitrarily complex non-linearity
with limited quantitative results. Only a few works provide characterizations
on the needed width or depth of GNN layers [li2020distance,
vignac2020building, loukas2020graph, zhang2021nested]. More quantitative
arguments on GNN analysis often depend on linear or Lipschitz continuous
assumptions to enable graph spectral analysis, such as feature oversmoothing
[oono2019graph, li2018deeper] and over-squashing [alon2020bottleneck,
topping2021understanding], the failure to process heterophilic graphs
[zhu2020beyond, yan2021two, chien2021adaptive] and the limited spectral
representation [balcilar2020analyzing, bianchi2021graph]. Some works also
study the generalization bounds [du2019graph, garg2020generalization,
liao2020pac] and the stability of GNNs [gama2019diffusion, gama2019stability,
levie2021transferability, gama2020stability]. However, the obtained results
may not reveal a direct comparison between non-linearity and linearity of the
model, and their analytic techniques avoid tackling the specific forms of non-
linear activations by using a Lipschitz continuous bound which is too loose in
our case.

Stochastic block models (SBM) and its contextual counterparts have been widely
used to study the node classification problems [abbe2017community,
abbe2018proof, massoulie2014community, bordenave2018nonbacktracking,
montanari2016semidefinite, binkiewicz2017covariate, deshpande2018contextual],
while these studies focus on the fundamental limits. Recently, (C)SBM and its
large-graph limitation also have been used to study the transferrability and
expressive power of GNN models [keriven2020convergence,
keriven2021universality, ruiz2020graphon] and GNNs on line graphs
[chen2020supervised], while these works did not compare non-linear and linear
propagation. CSBM has also been used to show the advantage of linear
convolution over no convolution for node classification [baranwal2021graph]. A
very recent result shows that attention-based propagation
[velivckovic2018graph] may be much worse than linear propagation given low-
quality node attributes under CSBM [fountoulakis2022graph]. Our results imply
that ReLU is the de facto optimal non-linearity instead of attention and may
at most marginally outperform the linear model when with low-quality node
attributes. Some previous works also use Bayesian inference to inspire GNN
architectures [jin2019graph, qu2019gmnn, kuck2020belief, jia2022unifying,
satorras2021neural, jia2021graph, qu2021neural, mehta2019stochastic], while
these works focus on empirical evaluation instead of theoretical analysis.

##  2 Preliminaries

In this section, we introduce preliminaries and notations for our later
discussion.

Maximum-a-posteriori (MAP) estimation. Suppose there are a set of finite
classes $\mathcal{C}$. A class label $Y\in\mathcal{C}$ is generated with
probability $\pi_{Y}$, where $\sum_{Y\in\mathcal{C}}\pi_{Y}=1$. Given $Y$, the
corresponding feature $X$ in the space $\mathcal{X}$ is generated from the
distribution $X\sim\mathbb{P}_{Y}$. A classifier is a decision
$f:\mathcal{X}\rightarrow\mathcal{C}$ and the _Bayesian mis-classification
error_ can be characterized as $\xi(f)=\sum_{Y\in\mathcal{C}}\pi_{Y}\int
1_{f(X)\neq Y}\mathbb{P}_{Y}(X)$, where and later $1_{S}$ indicates 1 if $S$
is true and 0 otherwise. The _MAP estimation_ of $Y$ given $X$ is the
classifier
$f^{*}(X)\triangleq\arg\max_{Y\in\mathcal{C}}\pi_{Y}\mathbb{P}_{Y}(X)$ that
can minimize $\xi(f)$ [theodoridis2015machine]. Later, we denote the minimal
Bayesian mis-classification error $\xi(f)$ as $\xi^{*}=\xi(f^{*})$.

Signal-to-Noise Ratio (SNR). Detection of a signal from the background
essentially corresponds to a binary classification problem. SNR is widely used
to measure the potential detection performance before specifying the
classifier [poor2013introduction]. In particular, if we have two equiprobable
classes $\mathcal{C}=\\{-1,1\\}$ and the features follows 1-d Gaussian
distributions $\mathbb{P}_{Y}=\mathcal{N}(\mu_{Y},\sigma^{2})$,
$Y\in\mathcal{C}$. The SNR $\rho$ defined as follows precisely characterizes
the minimal Bayesian mis-classification error.

$\displaystyle\textbf{SNR:}\quad\rho=\frac{\text{mean
difference}^{2}}{\text{variance}}=\frac{(\mu_{1}-\mu_{-1})^{2}}{\sigma^{2}}.$
(1)

In this case, the MAP estimation
$f^{*}(X)=2*1_{|X-\mu_{1}|\geq|X-\mu_{-1}|}-1$ and the minimal Bayesian mis-
classification error is $\Phi(-\sqrt{\rho}/2)$ where $\Phi$ denotes the
cumulative standard Gaussian distribution function. For more general cases
where the two classes are associated with sub-Guassian distributions
$\mathbb{P}_{Y}$, s.t.
$\mathbb{P}_{Y}(|X|>t)\in[c_{1}\exp(-c_{2}t^{2}),C_{1}\exp(-C_{2}t^{2})]$, for
some non-negative constants $c_{1},c_{2},C_{1},C_{2}$, a similar connection
between $\xi^{*}$ and $\rho$ can be shown by leveraging sharp sub-Gaussian
lower bounds [zhang2020non]. We will specify the connection to SNR in our case
in Sec. 4 and the SNR $\rho$ will be used as the main bridge to compare the
mis-classification errors of non-linear v.s. linear models.

Contextual Stochastic Block Model (CSBM). Random graph models have been widely
used to study the performance of algorithms on graphs [ding2021efficient,
li2019optimizing]. For node classification problems, CSBM is often used
[keriven2020convergence, keriven2021universality, ruiz2020graphon], as it well
combines the models of network structure and node attributes.

We study the case that nodes are in two equi-probable classes
$\mathcal{C}=\\{-1,1\\}$, where $\pi_{Y}=\frac{1}{2},Y\in\mathcal{C}$. Our
analysis can be generalized. An attributed network
$\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{X})$ is sampled from CSBM with
parameters $(n,p,q,\mathbb{P}_{1},\mathbb{P}_{-1})$ as follows. Suppose there
are $n$ nodes, $\mathcal{V}=[n]$. For each node $v$, the label
$Y_{v}\in\mathcal{C}$ is sampled from Rademacher distribution. Given $Y_{v}$,
the node attribute $X_{v}$ is sampled from $\mathbb{P}_{Y_{v}}$. For two nodes
$u,v$, if $Y_{u}=Y_{v}$, there is an edge $e_{uv}\in\mathcal{E}$ connecting
them with probability $p$. If $Y_{u}\neq Y_{v}$, there is an edge
$e_{uv}\in\mathcal{E}$ connecting them with probability $q$. All node
attributes $\mathbf{X}$ and edges $\mathcal{E}$ are independent given the node
labels $\mathbf{Y}=\\{Y_{v}|v\in\mathcal{V}\\}$.

Note that $p>q$ indicates the nodes with the same labels tend to be directly
connected, which corresponds to the _homophilic_ case, while $p<q$ corresponds
to the _heterophilic_ case.

The gap $|p-q|$, representing probabilities difference of a node connects to
nodes from the same class or the different class, reflects _structural
information_ and the gap between $\mathbb{P}_{1},\,\mathbb{P}_{-1}$ reflects
_attributed information_ , e.g., Jensen-Shannon distance
$\text{JS}(\mathbb{P}_{1},\mathbb{P}_{-1})$ that is well connected to Bayesian
mis-classification error [lin1991divergence]. Graph learning allows combining
these two types of information. In Sec. 4, we give more specific definitions
of these two types of information and their regime for our analysis.

##  3 Bayesian Inference and Nonlinearity in Graph Neural Networks

In the previous section, we discuss that given conditioned feature
distributions $X\sim\mathbb{P}_{Y}$, $Y\in\mathcal{C}$, the MAP estimation
$f^{*}(X)$ can minimize mis-classification error. For node classification in
an attributed network, the estimation of a node label should depend on not
only one’s own attributes but also its neighbors’. For example, in a
homophilic network, nodes with same labels tend to be directly connected.
Intuitively, using the averaged neighbor attributes may provide better
estimation of the label, which gives us graph convolution. In a heterophilic
network, nodes with different labels tend to be directed connected. So,
intuitively, checking the difference between one’s attributes and the
neighbors’ may provide better estimation. However, what could be the optimal
form to combine one’s own attribute with the neighbors’ attributes? We resort
to the MAP estimation. That is, given the attributes of a node
$v\in\mathcal{V}$ and its neighbors $\mathcal{N}_{v}$, we consider the MAP
estimation as follows.

$\displaystyle
f^{*}(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}})=\operatorname*{argmax}_{Y_{v}\in\mathcal{C}}\max_{Y_{u}\in\mathcal{C},\forall
u\in\mathcal{N}_{v}}\pi_{Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}}\mathbb{P}\left(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}},\mathcal{N}_{v}|Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}\right),$

where $\pi_{Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}}$ denotes their prior
distributions of node labels. Note that here we simplify the problem and
consider only 1-hop neighbors by following the setting [baranwal2021graph]. In
practice, most GNN models can only work on local networks due to the
scalability constraints [hamilton2017inductive, zeng2019graphsaint,
yin2022algorithm]. Even with the above simplification, the above MAP
estimation is generally intractable.

Therefore, we consider the CSBM with parameters
$(n,p,q,\mathbb{P}_{1},\mathbb{P}_{-1})$. In this case, the prior distribution
follows
$\pi_{Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}}=2^{-|\mathcal{N}_{v}|-1}$,
which is a constant given $\mathcal{N}_{v}$. The rest term follows

$\displaystyle\mathbb{P}\left(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}},\mathcal{N}_{v}|Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}\right)$
$\displaystyle=\mathbb{P}\left(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}}|Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}\right)\mathbb{P}\left(\mathcal{N}_{v}|Y_{v},\\{Y_{u}\\}_{u\in\mathcal{N}_{v}}\right)$
$\displaystyle=\prod_{u\in\mathcal{N}_{v}\cup\\{v\\}}\mathbb{P}_{Y_{u}}\left(X_{u}\right)\prod_{u\in\mathcal{N}_{v}}p^{(1+Y_{v}Y_{u})/2}q^{(1-Y_{v}Y_{u})/2}$
(2)

Therefore, the MAP estimation $f^{*}(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}})$
is to solve

$\displaystyle
f^{*}(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}})=\operatorname*{argmax}_{Y_{v}\in\mathcal{C}}\;\mathbb{P}_{Y_{v}}\left(X_{v}\right)\prod_{u\in\mathcal{N}_{v}}\max_{Y_{u}\in\mathcal{C}}\;\mathbb{P}_{Y_{u}}\left(X_{u}\right)p^{(1+Y_{v}Y_{u})/2}q^{(1-Y_{v}Y_{u})/2}$
(3)

This can be solved via the max-product algorithm [weiss2001optimality]. To
establish the connection to GNNs, we rewrite the RHS of Eq. 3 in the
logarithmic form and use the fact that $\mathcal{C}=\\{-1,1\\}$. And, we
achieve

$\displaystyle
f^{*}(X_{v},\\{X_{u}\\}_{u\in\mathcal{N}_{v}})=\text{sgn}\left(\log\frac{\mathbb{P}_{1}\left(X_{v}\right)}{\mathbb{P}_{-1}\left(X_{v}\right)}+\sum_{u\in\mathcal{N}_{v}}\mathcal{M}(X_{u},p,q)\right),\quad\text{where}$
$\displaystyle\mathcal{M}(X_{u},p,q)=\text{ReLU}\left(\log\frac{\mathbb{P}_{1}\left(X_{v}\right)}{\mathbb{P}_{-1}\left(X_{v}\right)}+\log\frac{p}{q}\right)-\text{ReLU}\left(\log\frac{\mathbb{P}_{1}\left(X_{v}\right)}{\mathbb{P}_{-1}\left(X_{v}\right)}+\log\frac{q}{p}\right)+\log\frac{q}{p}.$

We leave the derivation in Appendix LABEL:app:pro. Amazingly, activation ReLUs
in the message $\mathcal{M}$ well connect to the activations commonly-used in
GNN models, e.g., graph convolution networks [kipf2016semi]. Given the
optimality of the MAP estimation, we summarize this observation in Proposition
1.

######  Proposition 1 (Optimal Nonlinear Propagation).

Consider a network $\mathcal{G}\sim$
CSBM$(n,p,q,\mathbb{P}_{1},\mathbb{P}_{-1})$. To classify a node $v$, the
optimal nonlinear propagation (derived by the MAP estimation) given the
attributes of $v$ and its neighbors follows:

$\displaystyle\mathcal{P}_{v}=\psi\left(X_{v};\mathbb{P}_{1},\mathbb{P}_{-1}\right)+\sum_{u\in\mathcal{N}_{v}}\phi\left(\psi\left(X_{u};\mathbb{P}_{1},\mathbb{P}_{-1}\right);\log(p/q)\right)$
(4)

where
$\psi\left(a;\mathbb{P}_{1},\mathbb{P}_{-1}\right)=\log\frac{\mathbb{P}_{1}\left(a\right)}{\mathbb{P}_{-1}\left(a\right)}$
and
$\phi(a;\log\frac{p}{q})=\text{ReLU}(a+\log\frac{p}{q})-\text{ReLU}(a-\log\frac{p}{q})-\log\frac{p}{q}.$

Figure 2: Function $\phi(x;\log\frac{p}{q})$.

The optimal nonlinear propagation in Eq. (4) may contain two types of non-
linear functions: (1) $\psi$ is to measure the likelihood ratio between two
classes given the node attributes; (2) $\phi$ is to propagate the likelihood
ratios of the neighbors. ReLUs in $\phi$ avoid the overuse of the likelihood
ratios from neighbors, as $\phi$ essentially provides a bounded function (See
Fig. 2). One observation of the direct benefit of this non-linear propagation
is as follows.

######  Remark 1.

When there is no structural information, i.e., $p=q$, $\phi(x;0)=0$, $\forall
x\in\mathbb{R}$, the propagation is deactivated, which avoids potential
contamination from the attributes of the neighbors.

In the equiprobable case, the MAP estimation also gives the maximum likelihood
estimation (MLE) of $Y$ if we view the labels as the fixed parameters. When
the classes are unbalanced $\pi_{Y}\neq\frac{1}{2}$, similar results can be
obtained while additional terms $\log\frac{\pi_{1}}{\pi_{-1}}$ may appear as
bias in Eq. (4). Later, our analysis focuses on the equiprobable case while
empirical results in Sec. 5 show more general cases.

Moreover, if one is to infer the posterior distribution of $Y$, one may
replace the max-product algorithm to solve Eq. (3) with the sum-product
algorithm [pearl1982reverend]. Then, the obtained non-linearity in $\phi$ will
turn into Tanh functions. As ReLUs are more used in practical GNNs, we focus
on the case with ReLUs.

Discussion on the Non-linearity. Next, we discuss more insights into the non-
linearity of $\psi$ and $\phi$.

The function $\psi$ essentially corresponds to _a node-attribute
transformation_ , which depends on the distributions $\mathbb{P}_{\pm 1}$. As
these distributions are unknown in practice, a NN model to learn $\psi$ is
suggested, such as the one in the model APPNP [klicpera2018predict] and GPR-
GNN [chien2021adaptive]. Due to the expressivity of NNs
[cybenko1989approximation, hornik1989multilayer], a wide range of $\psi$ can
be modeled. One insightful example is that when $\mathbb{P}_{\pm 1}$ are
Laplace distributions, $\psi$ is a bounded function (same as $\phi$) to
control the heavy-tailed attributes generated by Laplace distributions.

######  Example 1 (Laplace Assumption).

When node attributes follow $m$-dimensional independent Laplace distribution,
i.e., $\mathbb{P}_{Y_{v}}=\frac{1}{(2b)^{m}}\exp(-\|X_{v}-Y_{v}\mu\|_{1}/b)$
for $\mu\in\mathbb{R}^{m}$, $b>0$ and $Y_{v}\in\\{-1,1\\}$. According to Eq.
(4), the function $\psi(\cdot;\mathbb{P}_{1},\mathbb{P}_{-1})$ can be
specified as

$\displaystyle\psi_{\text{lap}}(X_{v};\mathbb{P}_{1},\mathbb{P}_{-1})=\mathbbm{1}^{T}\phi(X_{v};2\mu/b),\;\text{where}\;\phi\;\text{as
defined in Eq.~{}\eqref{eq:prop} works in an entry-wise manner.}$

As node-attribute distributions may vary a lot, $\psi$ is better to be modeled
via a general NN in practice. More interesting findings may come from $\phi$
in Eq. (4) as it has a fixed form and well matches the most commonly-used GNN
architecture. Specifically, besides the extreme case stated in Remark 1, we
are to investigate _how_ non-linearity induced by the ReLUs in $\phi$ may
benefit the model. We expect the findings to provide the explanations to some
previous empirical experiences on using GNNs.

To simplify our discussion, when analyzing $\phi$, we focus on the case with a
linear node-attribute transformation $\psi=\psi_{\text{Gau}}$ in Eq. (6) by
assuming label-dependent Gaussian node attributes. This follows the
assumptions in previous studies [baranwal2021graph, fountoulakis2022graph]. In
fact, there are a class of distributions named natural exponential family
(NEF) [morris1982natural] which if the node attributes satisfy, the induced
$\phi$ is linear. We conjecture that our later results in Sec. 4 are applied
to the general NEF since the only difference is the bias term by comparing Eq.
(5) and Eq. (6).

######  Example 2 (Natural Exponential Family Assumption).

When node features follow $m$-dimensional natural exponential family
distributions
$\mathbb{P}_{Y_{v}}(X)=h(X_{v})\cdot\exp(\theta_{Y_{v}}^{T}X_{v}-M(\theta_{Y_{v}}))$
for $\theta_{Y_{v}}\in\mathbb{R}^{m}$ and $Y_{v}\in\\{-1,1\\}$ where
$M(\theta_{Y_{v}})$ is a parameter function. The function
$\psi(\cdot;\mathbb{P}_{1},\mathbb{P}_{-1})$ is specified as:

$\displaystyle\psi_{\text{nef}}(X_{v};\theta_{1},\theta_{-1})=(\theta_{1}-\theta_{-1})^{T}X_{v}-(M(\theta_{1})-M(\theta_{-1})).$
(5)

In particular, when $\mathbb{P}_{1}=\mathcal{N}(\mu,I/m)$,
$\mathbb{P}_{-1}=\mathcal{N}(\nu,I/m)$ for $\mu,\nu\in\mathbb{R}^{m}$,

$\displaystyle\psi_{\text{Gau}}(X_{v};\mu,\nu)=m\left[(\mu-\nu)^{T}X_{v}-(\|\mu\|_{2}^{2}-\|\nu\|_{2}^{2})/2\right].$
(6)

More generally, our optimal nonlinear propagation Eq. (4) can be well
generalized to other settings as long as the model satisfies edge-independent
assumption, where edges random variables are mutually independent conditioned
on the labels of nodes. When this assumption is satisfied, the MAP estimation
will result in graph convolution with ReLU activation.

We summarize our main theoretical findings regarding the nonlinearity of
$\phi$ in the next section.

##  4 Main Results on ReLU-based Nonlinear Propagation

In this section, we summarize our analytical results on $\phi$ in the optimal
nonlinear propagation (Eq. (4)). Our study assumes an attributed network
generated from CSBM$(n,p,q,\mathcal{N}(\mu,I/m),\mathcal{N}(\nu,I/m))$ where
$\mu,\nu\in\mathbb{R}^{m}$. We use CSBM-G$(n,p,q,\mu,\nu)$ later to denote
this model for simplicity. We are interested in the asymptotic behavior when
$n\rightarrow\infty$. Note that all parameters $\mu,\nu,p,q,m$ may implicitly
depend on $n$. We are to compare the non-linear propagation model
$\mathcal{P}_{v}$ suggested by Eq. (4) where $\psi=\psi_{\text{Gau}}$ with the
following linear counterpart $\mathcal{P}_{v}^{l}$.

$\displaystyle\textbf{Baseline linear model:}\quad\mathcal{P}_{v}^{l}(w)$
$\displaystyle=\psi_{\text{Gau}}(X_{v};\mu,\nu)+w\sum_{u\in\mathcal{N}_{v}}\psi_{\text{Gau}}(X_{u};\mu,\nu),\;\text{for
all $v\in\mathcal{V}$.}$ (7)

where $w\in\mathbb{R}$ is an extra parameter to be tuned. Note that this
linear model can be claimed as an optimal linear model up-to a choice of $w$
because the distributions of both the center node attribute $X_{v}$ and the
linear aggregation from the neighbors $\sum_{u\in\mathcal{N}_{v}}X_{u}$ are
Gaussian and symmetric w.r.t. the hyperplane
$\\{Z\in\mathbb{R}^{m}|(\mu-\nu)^{T}Z=(\|\mu\|_{2}^{2}-\|\nu\|_{2}^{2})/2\\}$
for the two classes. We are to compare their classification errors
$\xi^{r}=\xi(\text{sgn}(\mathcal{P}_{v}))$ and
$\xi^{l}(w)=\xi(\text{sgn}(\mathcal{P}_{v}^{l}(w)))$. By following
[baranwal2021graph], we also discuss separability of all nodes in the network,
i.e., $\mathbb{P}(\forall v\in\mathcal{V},\mathcal{P}_{v}\cdot Y_{v}>0)$ in
Theorem 1.

To begin with, we introduce several quantities for the convenience of further
statements. The SNRs

$\rho_{r}=\frac{\left(\mathbb{E}[\mathcal{P}_{v}|Y_{v}=1]-\mathbb{E}[\mathcal{P}_{v}|Y_{v}=-1]\right)^{2}}{\text{var}(\mathcal{P}_{v}|Y_{v}=1)},\rho_{l}(w)=\frac{\left(\mathbb{E}[\mathcal{P}_{v}^{l}(w)|Y_{v}=1]-\mathbb{E}[\mathcal{P}_{v}^{l}(w)|Y_{v}=-1]\right)^{2}}{\text{var}(\mathcal{P}_{v}^{l}(w)|Y_{v}=1)}$

are important quantities to later characterize different types of propagation.
Also, we characterize structural information by
$\mathcal{S}(p,q)=(p-q)^{2}/(p+q)$ and attributed information by
$\sqrt{m}\|\mu-\nu\|_{2}$.

######  Assumption 1 (Moderate Structural Information).

$\mathcal{S}(p,q)=\omega_{n}(\frac{(\log n)^{2}}{n})$ and
$\frac{\mathcal{S}(p,q)}{|p-q|}\nrightarrow 1$.

######  Assumption 2 (Bounded Attributed Information).

$\sqrt{m}\|\mu-\nu\|_{2}=o_{n}(\log n)$.

Assumption 1 states that structural information should be neither too weak nor
too strong. $\mathcal{S}(p,q)=\omega_{n}(\frac{(\log n)^{2}}{n})$ excludes the
extremely weak case discussed in Remark 1. Moreover, the graph structure
should not be too sparse, so the aggregated information from neighbors
dominates the propagation. $\frac{\mathcal{S}(p,q)}{|p-q|}\nrightarrow 1$
means neither $p=\omega_{n}(q)$ nor $q=\omega_{n}(p)$, which avoids extremely
strong structural information. This assumption is more general than some
concurrent works on CSBM-based GNN analysis [baranwal2021graph,
fountoulakis2022graph] as we include the cases with less structural
information $|p-q|=o_{n}(p+q)$ and with heterophily $p<q$. Assumption 2 is to
avoid too strong attributed information: when
$\sqrt{m}\|\mu-\nu\|_{2}=\Omega_{n}(\log n)$, all nodes in CSBM can be
accurately classified in the asymptotic sense without structural information,
i.e. $\mathbb{P}(\forall
v\in\mathcal{V},\,\psi_{\text{Gau}}(X_{v};\mu,\nu)\cdot Y_{v}>0)=1-o_{n}(1)$.
Now, we present our first lemma which links the mis-classification errors
$\xi^{r},\xi^{l}$ to with the SNRs $\rho_{r},\rho_{l}$:

######  Lemma 1.

Suppose $(p,q)$ satisfies Assumption 1, for any $\mathcal{G}\sim$
CSBM-G$(n,p,q,\mu,\nu)$,

$\displaystyle\xi^{r}\in[C_{1}\exp(-C_{2}\rho_{r}/2),\exp(-\rho_{r}/2)],\quad\xi^{l}(w)\to\exp(-\rho_{l}(w)(1+o_{n}(1))/2)$
(8)

where $C_{2}$ is asymptotically a constant, and the notation $a(n)\to b(n)$
denotes $a(n)/b(n)\to 1$. Lemma 1 claims that the classification errors under
both nonlinear and linear model can be controlled by their SNRs. By leveraging
Lemma 1, we can further illustrate the separability of all nodes in the
network, which is presented in the following theorem.

######  Theorem 1 (Separability).

Suppose that $(p,q)$ satisfies Assumption 1, for $\mathcal{G}\sim$
CSBM-G$(n,p,q,\mu,\nu)$, if $\sqrt{m}\|\mu-\nu\|_{2}=\omega_{n}(\sqrt{\log
n/S(p,q)n})$, then

$\displaystyle\mathbb{P}(\forall v\in\mathcal{V},\mathcal{P}_{v}\cdot
Y_{v}>0)=1-\mathcal{O}_{n}(n\exp(-\rho_{r}/2))=1-o_{n}(1),$ (9)
$\displaystyle\mathbb{P}(\forall v\in\mathcal{V},\mathcal{P}_{v}^{l}(w)\cdot
Y_{v}>0)=1-\mathcal{O}_{n}(n\exp(-\rho_{l}(w)/2))=1-o_{n}(1).$ (10)

Here, assume $|w|>c$ for some positive constant $c$ and
$\text{sgn}(w)=\text{sgn}(p-q)$ in the linear model.

Theorem 1 applies to both homophilic ($p>q$) and heterophilic $(p<q)$
scenarios. Even for just the linear case, compared to [baranwal2021graph]
which needs $\sqrt{m}\|\mu-\nu\|_{2}=\omega_{n}(\log n/\sqrt{S(p,q)n})$ to
achieve separability, we have $\sqrt{\log n}$ improvement due to a tight
analysis.

As shown in Lemma 1 and Theorem 1, the errors are mainly determined by SNRs.
Large SNR implies a fast decay rate of the errors of a single node and the
entire graph, which motivates us to further explore SNRs to illustrate a
comparison between non-linear and linear models. We consider comparing with
the optimal linear model, i.e., $\rho_{l}^{*}=\rho_{l}(w^{*})$, where
$w^{*}=\arg\min_{w\in\mathbb{R}}\xi^{l}(w)$.

######  Theorem 2.

Suppose that $(p,q)$ satisfies Assumption 1, for $\mathcal{G}\sim$
CSBM-G$(n,p,q,\mu,\nu)$, under the separable condition in Theorem 1
$\sqrt{m}\|\mu-\nu\|_{2}=\omega_{n}(\sqrt{\log n/S(p,q)n})$, we further have

  * •

I. Limited Attributed Information: When
$\sqrt{m}\|\mu-\nu\|_{2}=\mathcal{O}_{n}(1)$,

$\displaystyle\rho_{r}=\Theta_{n}(\rho_{l}^{*}),$ (11)

Further, if $\sqrt{m}\|\mu-\nu\|_{2}=o_{n}(|\log(p/q)|)$,
$\rho_{r}/\rho_{l}^{*}\to 1$;

  * •

II. Sufficient Attributed Information: When
$\sqrt{m}\|\mu-\nu\|_{2}=\omega_{n}(1)$ and satisfies Assumption 2,

$\displaystyle\rho_{r}=$
$\displaystyle\omega_{n}(\min\\{\exp(m\|\mu-\nu\|_{2}^{2}/3),nS(p,q)m^{-1}\|\mu-\nu\|_{2}^{-2}\\}\cdot\rho_{l}^{*})=\omega_{n}(\rho_{l}^{*}).$
(12)

Theorem 2 also works both homophilic ($p>q$) and heterophilic $(p<q)$
scenarios. Theorem 2 implies that when attributed information is limited,
nonlinear propagation behaves similar to the linear model as their SNRs are in
the same order. Particularly, when attributed information is very limited,
$\sqrt{m}\|\mu-\nu\|_{2}=o_{n}(|\log(p/q)|)$, the SNRs of two models are
asymptotically the same. In the regime of sufficient attributed information,
nonlinear propagation brings order-level superiority compared with the linear
model. The intuition is that in this regime, when the attributes are very
informative, the bounds of $\phi$ in Eq. (4) help with avoiding overconfidence
given by the node attributes. The coefficient before $\rho_{l}^{*}$ in Eq.
(12) shows the trade-off between structural information and attributed
information on controlling the superiority of nonlinear propagation.

Next, we analyze whether nonlinearity makes model more transferable or not
when there often exists a distribution shift between the training and testing
datasets, which is also practically useful.

We consider the following setting. We assume using a large enough network
generated by CSBM-G$(n,p,q,\mu,\nu)$ for training so that the optimal
parameters as in $\mathcal{P}_{v}$ and $\mathcal{P}_{v}^{l}(w^{*})$ for this
CSBM-G have been learnt. We consider their mis-classification errors over
another CSBM-G with parameters
$(n,p^{\prime},q^{\prime},\mu^{\prime},\nu^{\prime})$. We keep the amounts of
attributed information and structural information unchanged by setting
$p=p^{\prime}$, $q=q^{\prime}$,
$\mu^{\prime}=(\mu+\nu)/2+\bm{R}(\mu-\nu)/2,\nu^{\prime}=\mu+\nu/2+\bm{R}(\nu-\mu)/2$
for a rotation matrix $\bm{R}$ close to $I$. Let $\Delta\xi^{r}$ and
$\Delta\xi{{}^{l}}(w^{*})$ denote the increase of mis-classification errors of
models $\mathcal{P}_{v}$ and $\mathcal{P}_{v}^{l}(w^{*})$, respectively, due
to such a distribution shift. We may achieve the following results.

######  Theorem 3 (Transferability).

Suppose that $(p,q)$ satisfies Assumption 1, for $\mathcal{G}^{\prime}\sim$
CSBM-G$(n,p^{\prime},q^{\prime},\mu^{\prime},\nu^{\prime})$, under the linear
separable condition
$\sqrt{m}\|\mu^{\prime}-\nu^{\prime}\|_{2}=\omega_{n}(\sqrt{\log
n/S(p^{\prime},q^{\prime})n})$. Suppose $\mathcal{P}_{v}$ and
$\mathcal{P}_{v}^{l}(w^{*})$ have learnt parameters from $\mathcal{G}\sim$
CSBM-G$(n,p,q,\mu,\nu)$ where the parameters of two CSBM-Gs follow the
relation described above. Then, we have

  * •

I. Limited Attributed Information: When
$\sqrt{m}\|\mu-\nu\|_{2}=o_{n}(|\log(p/q)|)$,
$\Delta\xi^{r}/\Delta\xi^{l}(w^{*})\to 1$.

  * •

II. Sufficient Attributed Information: When
$\sqrt{m}\|\mu-\nu\|_{2}=\omega_{n}(1)$ and satisfies Assumption 2,
$\Delta\xi^{r}/\Delta\xi{{}^{l}}(w^{*})\to 0$.

Similar to Theorem 2, when attributed information is very limited,
nonlinearity will not bring any benefit, while in the regime with informative
attributes, nonlinearity increases model transferability. We leave the
intermediate regime
$\sqrt{m}\|\mu-\nu\|_{2}\in[\Omega_{n}(|\log(p/q)|),O_{n}(1)]$ for future
study.

##  5 Experiments

In this section, we verify our theoretical results based on synthetic and real
datasets. In all experiments, we fix $w$ in the linear model (Eq. (7)) as
$w=1$ for the homophilic case $(p>q)$ and $w=-1$ for the heterophilic case
$(p<q)$. Experiments on other $w$’s can be found in Appendix LABEL:app:exp_w,
which does not change the achieved conclusion. This is because when the node
number $n$ is large, for a constant $w$, the neighbor information will
dominate the results. Later, we use
$\mathcal{P}_{v}^{l}=\mathcal{P}_{v}^{l}(w)$ for simplicity.

###  5.1 Asymptotic Experiments - Model Accuracy & Transferability Study

Our first experiments focus on evaluating the asymptotic
$(n\rightarrow\infty)$ classification performance of nonlinear and linear
models. Given a CSBM-G, we generate 5 graphs and compute the average accuracy
results (#correctly classified nodes / #total nodes). We compare the nonlinear
v.s. linear models under three different CSBM-G settings. Fig. 3 shows the
results.

Figure 3: Classification Performance on Nonlinear Models v.s. Linear Models
($\mathcal{P}_{v}$ v.s. $\mathcal{P}_{v}^{l}$). LEFT: Homophily + Limited
Attr. Info.; MIDDLE: Homophily + Suff. Attr. Info.; RIGHT: Heterophily + Fixed
Attr. Info.. $m=10$ and other parameters are listed in the figures.

All three cases satisfy the separability condition in Theorem 1, so, as $n$
increases, the accuracy progressively increases to 1. Our results also match
well with the implications provided by Theorem 2. In the regime with limited
attribute information (Fig. 1 LEFT) where $\rho_{r}=\Theta_{n}(\rho_{l}^{*})$
as proved, the nonlinear model and the linear model behave almost the same
(performance gap $<0.15\%$ for $n\geq 10^{5}$). In the regime with sufficient
attribute information (Fig. 1 MIDDLE) where
$\rho_{r}=\omega_{n}(\rho_{l}^{*})$ as proved, we may observe that the
nonlinear model can significantly outperform the linear model as
$n\rightarrow\infty$. Fig. 1 RIGHT is to show the heterophilic graph case
($p<q$). If we switch the values of $p,q$ (and also change the models
correspondingly), we obtain the exactly same figure up to some experimental
randomness (see Appendix LABEL:app:heter_to_homo). Also, Fig. 1 RIGHT
considers a boundary case of sufficient attributed information, i.e.,
$\sqrt{m}\|\mu-\nu\|_{2}=\Theta_{n}(1)$. We observe that Theorem 2 still well
describes the asymptotic performance when $n\rightarrow\infty$.

We further study the transferability for the non-linear model and the linear
model. We follow the setting in Theorem 3 by rotating
$\mu,\nu\rightarrow\mu^{\prime},\nu^{\prime}$. Fig. LABEL:fig.perturb shows
the result and well matches Theorem 3. In the regime of limited attributed
information, the two models have the almost same transferability, i.e., the
perturbation error ratio is close to 1. In contrast, with sufficient
attributed information, the non-linear model is more transferrable than the
linear counterpart as the ratio is smaller than 1.

2.5cm(5.6cm,3.5cm)

$p=\frac{2}{\sqrt{n}}$

$q=\frac{1}{\sqrt{n}}$

$\|\mu-\nu\|_{2}=\frac{0.3\log^{2}n}{\sqrt{n}}$ 2.5cm(10.2cm,3.5cm)