GRAPH NEURAL NETWORKS EXPONENTIALLY LOSE EXPRESSIVE POWER FOR NODE CLASSIFICATION
Year: 2020

Abstract:
Graph Neural Networks (graph NNs) are a promising deep learning approach for analyzing graph-structured data. However, it is known that they do not improve (or sometimes worsen) their predictive performance as we pile up many layers and add non-lineality. To tackle this problem, we investigate the expressive power of graph NNs via their asymptotic behaviors as the layer size tends to infinity. Our strategy is to generalize the forward propagation of a Graph Convolutional Network (GCN), which is a popular graph NN variant, as a specific dynamical system. In the case of a GCN, we show that when its weights satisfy the conditions determined by the spectra of the (augmented) normalized Laplacian, its output exponentially approaches the set of signals that carry information of the connected components and node degrees only for distinguishing nodes. Our theory enables us to relate the expressive power of GCNs with the topological information of the underlying graphs inherent in the graph spectra. To demonstrate this, we characterize the asymptotic behavior of GCNs on the Erdős – Rényi graph. We show that when the Erdős – Rényi graph is sufficiently dense and large, a broad range of GCNs on it suffers from the “information loss” in the limit of infinite layers with high probability. Based on the theory, we provide a principled guideline for weight normalization of graph NNs. We experimentally confirm that the proposed weight scaling enhances the predictive performance of GCNs in real data1.

heading: 1 INTRODUCTION
text: Motivated by the success of Deep Learning (DL), several attempts have been made to apply DL models to non-Euclidean data, particularly, graph-structured data such as chemical compounds, social networks, and polygons. Recently, Graph Neural Networks (graph NNs) (Duvenaud et al., 2015; Li et al., 2016; Gilmer et al., 2017; Hamilton et al., 2017; Kipf & Welling, 2017; Nguyen et al., 2017; Schlichtkrull et al., 2018; Battaglia et al., 2018; Xu et al., 2019; Wu et al., 2019a) have emerged as a promising approach. However, despite their practical popularity, theoretical research of graph NNs has not been explored extensively.
The characterization of DL model expressive power, i.e., to identify what function classes DL models can (approximately) represent, is a fundamental question in theoretical research of DL. Many studies have been conducted for Fully Connected Neural Networks (FNNs) (Cybenko, 1989; Hornik, 1991; Hornik et al., 1989; Barron, 1993; Mhaskar, 1993; Sonoda & Murata, 2017; Yarotsky, 2017) and Convolutional Neural Networks (CNNs) (Petersen & Voigtlaender, 2018; Zhou, 2018; Oono & Suzuki, 2019). For such models, we have theoretical and empirical justification that deep and nonlinear architectures can enhance representation power (Telgarsky, 2016; Chen et al., 2018b; Zhou & Feng, 2018). However, for graph NNs, several papers have reported that node representations go indistinguishable (known as over-smoothing) and prediction performances severely degrade when we stack many layers (Kipf & Welling, 2017; Wu et al., 2019b; Li et al., 2018). Besides, Wu et al. (2019a) reported that graph NNs achieved comparable performance even if they removed intermediate non-linear functions. These studies posed a question about the current architecture and made us aware of the need for the theoretical analysis of the graph NN expressive power.
1Code is available at https://github.com/delta2323/gnn-asymptotics.
In this paper, we investigate the expressive power of graph NNs by analyzing their asymptotic behaviors as the layer size goes to infinity. Our theory gives new theoretical conditions under which neither layer stacking nor non-linearity contributes to improving expressive power. We consider a specific dynamics that includes a transition defining a Markov process and the forward propagation of a Graph Convolutional Network (GCN) (Kipf & Welling, 2017), which is one of the most popular graph NN variants, as special cases. We prove that under certain conditions, the dynamics exponentially approaches a subspace that is invariant under the dynamics. In the case of GCN, the invariant space is a set of signals that correspond to the lowest frequency of graph spectra and that have “no information” other than connected components and node degrees for a node classification task whose goal is to predict the nodes’ properties in a graph. The rate of the distance between the output and the invariant space is O((sλ)L) where s is the maximum singular values of weights, λ is typically a quantity determined by the spectra of the (augmented) normalized Laplacian, and L is the layer size. See Sections 3.3 (general case) and 4 (GCN case) for precise statements.
We can interpret our theorem as the generalization of the well-known property that if a finite and discrete Markov process is irreducible and aperiodic, it exponentially converges to a unique equilibrium and the eigenvalues of its transition matrix determine the convergence rate (see, e.g., Chung & Graham (1997)). Different from the Markov process case, which is linear, the existence of intermediate non-linear functions complicates the analysis. We overcame this problem by leveraging the combination of the ReLU activation function (Krizhevsky et al., 2012) and the positivity of eigenvectors of the Laplacian associated with the smallest positive eigenvalues.
Our theory enables us to investigate asymptotic behaviors of graph NNs via the spectral distribution of the underlying graphs. To demonstrate this, we take GCNs defined on the Erdős – Rényi graph GN,p, which has N nodes and each edge appears independently with probability p, for an example. We prove that if logNpN = o(1) as a function of N , any GCN whose weights have maximum singular
values at most C √
Np log(N/ε) approaches the “information-less” invariant space with probability at
least 1− ε, where C is a universal constant. Intuitively, if the graph on which we define graph NNs is sufficiently dense, graph-convolution operations mix signals on nodes fast and hence the feature maps lose information for distinguishing nodes quickly.
Our contributions are as follows:
• We relate asymptotic behaviors of graph NNs with the topological information of underlying graphs via the spectral distribution of the (augmented) normalized Laplacian.
• We prove that if the weights of a GCN satisfy conditions determined by the graph spectra, the output of the GCN carries no information other than the node degrees and connected components for discriminating nodes when the layer size goes to infinity (Theorems 1, 2).
• We apply our theory to Erdős – Rényi graphs as an example and show that when the graph is sufficiently dense and large, many GCNs suffer from the information loss (Theorem 3).
• We propose a principled guideline for weight normalization of graph NNs and empirically confirm it using real data.

heading: 2 RELATED WORK
text: MPNN-type Graph NNs. Since many graph NN variants have been proposed, there are several unified formulations of graph NNs (Gilmer et al., 2017; Battaglia et al., 2018). Our approach is the closest to the formulation of Message Passing Neural Network (MPNN) (Gilmer et al., 2017), which unified graph NNs in terms of the update and readout operations. Many graph NNs fall into this formulation such as Duvenaud et al. (2015), Li et al. (2016), and Veličković et al. (2018). Among others, GCN (Kipf & Welling, 2017) is an important application of our theory because it is one of the most widely used graph NNs. In addition, GCNs are interesting from a theoretical research perspective because, in addition to an MPNN-type graph NN, we can interpret GCNs as a simplification of spectral-type graph NNs (Henaff et al., 2015; Defferrard et al., 2016), that make use of the graph Laplacian.
Our approach, which considers the asymptotic behaviors graph NNs as the layer size goes to infinity, is similar to Scarselli et al. (2009), one of the earliest works about graph NNs. They obtained node
representations by iterating message passing between nodes until convergence. Their formulation is general in that we can use any local aggregation operation as long as it is a contraction map. Our theory differs from theirs in that we proved that the output of a graph NN approaches a certain space even if the local aggregation function is not necessarily a contraction map.
Expressive Power of Graph NNs. Several studies have focused on theoretical analysis and the improvement of graph NN expressive power. For example, Xu et al. (2019) proved that graph NNs are no more powerful than the Weisfeiler – Lehman (WL) isomorphism test (Weisfeiler & A.A., 1968) and proposed a Graph Isomorphism Network (GIN), that is approximately as powerful as the WL test. Although they experimentally showed that GIN has improved accuracy in supervised learning tasks, their analysis was restricted to the graph isomorphism problem. Xu et al. (2018) analyzed the non-asymptotic properties of GCNs through the lens of random walk theory. They proved the limitations of GCNs in expander-like graphs and proposed a Jumping Knowledge Network (JKNet) to address the issue. To handle the non-linearity, they linearized networks by a randomization assumption (Choromanska et al., 2015). We take a different strategy and make use of the interpretation of ReLU as a projection onto a cone. Recently, NT & Maehara (2019) showed that a GCN approximately works as a low-pass filter plus an MLP in a certain setting. Although they analyzed finite-depth GCNs, our theory has similar spirits with theirs because our “information-less” space corresponds to the lowest frequency of a graph Laplacian. Another point is that they imposed assumptions that input signals consist of low-frequent true signals and high-frequent noise, whereas we need not such an assumption.
Role of Deep and Non-linear Structures. For ordinal DL models such as FNNs and CNNs, we have both theoretical and empirical justification of deep and non-linear architectures for enhancing of the expressive power (e.g., Telgarsky (2016); Petersen & Voigtlaender (2018); Oono & Suzuki (2019)). In contrast, several studies have witnessed severe performance degradation when stacking many layers on graph NNs (Kipf & Welling, 2017; Wu et al., 2019b). Li et al. (2018) reported that feature vectors on nodes in a graph go indistinguishable as we increase layers in several tasks. They named this phenomenon over-smoothing. Regarding non-linearity, Wu et al. (2019a) empirically showed that graph NNs achieve comparable performance even if we omit intermediate non-linearity. These observations gave us questions about the current models of deep graph NNs in terms of their expressive power. Several studies gave theoretical explanations of the over-smoothing phenomena for linear GNNs (Li et al., 2018; Zhang, 2019; Zhao & Akoglu, 2020). We can think of our theory as an extension of their results to non-linear GNNs.

heading: 3 PROBLEM SETTING AND MAIN RESULT

heading: 3.1 NOTATION
text: Let N+ be the set of positive integers. For N ∈ N+, we denote [N ] := {1, . . . , N}. For a vector v ∈ RN , we write v ≥ 0 if and only if vn ≥ 0 for all n ∈ [N ]. Similarly, for a matrix X ∈ RN×C , we write X ≥ 0 if and only if Xnc ≥ 0 for all n ∈ [N ] and c ∈ [C]. We say such a vector and matrix is non-negative. 〈·, ·〉 denotes the inner product of vectors or matrices, depending on the context: 〈u, v〉 := u>v for u, v ∈ RN and 〈X,Y 〉 := tr(XTY ) for X,Y ∈ RN×C . 1P equals to 1 if the proposition P is true else 0. For vectors v ∈ RN and w ∈ RC , v ⊗ w ∈ RN×C denotes the Kronecker product of v and w defined by (v⊗w)nc := vnwc. For X ∈ RN×C , ‖X‖F := 〈X,X〉 1 2 denotes the Frobenius norm of X . For a vector v ∈ RN , diag(v) := (vnδnm)n,m∈[N ] ∈ RN×N denotes the diagonalization of v. IN ∈ RN×N denotes the identity matrix of size N . For a linear operator P : RN → RM and a subset V ⊂ RN , we denote the restriction of P to V by P |V .

heading: 3.2 DYNAMICAL SYSTEM
text: Although we are mainly interested in GCNs, we develop our theory more generally using dynamical systems. We will specialize to the GCNs in Section 4.
For N,C,Hl ∈ N+ (l ∈ N+), let P ∈ RN×N be a symmetric matrix and Wlh ∈ RC×C for l ∈ N+ and h ∈ [Hl]. We define fl : RN×C → RN×C by fl(X) := MLPl(PX). Here, MLPl : RN×C → RN×C is the l-th multi-layer perceptron common to all nodes (Xu et al., 2019) and is defined by MLPl(X) := σ(· · ·σ(σ(X)Wl1)Wl2 · · ·WlHl), where σ : RN×C → RN×C
is an element-wise ReLU function (Krizhevsky et al., 2012) defined by σ(X)nc := max(Xnc, 0) for n ∈ [N ], c ∈ [C]. We consider the dynamics X(l+1) := fl(X(l)) with some initial value X(0) ∈ RN×C . We are interested in the asymptotic behavior of X(l) as l→∞. For M ≤ N , let U be a M -dimensional subspace of RN . We assume that U and P satisfy the following properties that generalize the situation where U is the eigenspace associated with the smallest eigenvalue of a (normalized) graph Laplacian ∆ (that is, zero) and P is a polynomial of ∆. Assumption 1. U has an orthonormal basis (em)m∈[M ] that consists of non-negative vectors. Assumption 2. U is invariant under P , i.e., if u ∈ U , then Pu ∈ U .
We endow RN with the ordinal inner product and denote the orthogonal complement of U by U⊥ := {u ∈ RN | 〈u, v〉 = 0,∀v ∈ U}. By the symmetry of P , we can show that U⊥ is invariant under P , too (Appendix E.1, Proposition 2). Therefore, we can regard P as a linear mapping P |U⊥ : U⊥ → U⊥. We denote the operator norm of P |U⊥ by λ. When U is the eigenspace associated with the smallest eigenvalue of ∆ and P is g(∆) where g is a polynomial, then, λ corresponds to λ = supµ |g(µ)| where sup ranges over all eigenvalues except the smallest one.

heading: 3.3 MAIN RESULT
text: We define the subspaceM of RN×C byM := U ⊗ RC = { ∑M m=1 em ⊗ wm | wm ∈ RC} where (em)m∈[M ] is the orthonormal basis of U appeared in Assumption 1. For X ∈ RN×C , we denote the distance betweenX andM by dM(X) := inf{‖X−Y ‖F | Y ∈M}. We denote the maximum singular value of Wlh by slh and set sl := ∏Hl h=1 slh. With these preparations, we introduce the main theorem of the paper. Theorem 1. Under Assumptions 1 and 2, we have dM(fl(X)) ≤ slλdM(X) for any X ∈ RN×C .
The proof key is that the non-linear operation σ decreases the distance dM, that is, dM(σ(X)) ≤ dM(X). We use the non-negativity of em to prove this claim. See Appendix A for the complete proof. We also discuss the strictness of Theorem 1 in Appendix E.3.
By setting dM(X) = 0, this theorem implies that M is invariant under fl. In addition, if the maximum value of singular values are small, X(l) asymptotically approaches M in the sense of Johnson (1973) for any initial value X(0). That is, the followings hold under Assumptions 1 and 2. Corollary 1. M is invariant under fl for any l ∈ N+, that is, ifX ∈M, then we have fl(X) ∈M. Corollary 2. Let s := supl∈N+ sl. We have dM(X
(l)) = O((sλ)l). In particular, if sλ < 1, then Xl exponentially approachesM as l→∞ for any initial value X(0).
Suppose the operator norm of P |U : U → U is no larger than λ, then, under the assumption of sλ < 1, X(l) converges to 0, the trivial fixed point (see Appendix E.2, Proposition 3). Therefore, we are mainly interested in the case where the operator norm of P |U is strictly larger than λ (see Proposition 1). Finally, we restate Theorem 1 specialized to the situation where U is the direct sum of eigenspaces associated with the largest M eigenvalues of P . Note that the eigenvalues of P is real since P is symmetric. Corollary 3. Let λ1 ≤ · · · ≤ λN be the eigenvalue of P , sorted in ascending order. Suppose the multiplicity of the largest eigenvalue λN is M(≤ N), i.e., λN−M < λN−M+1 = · · · = λN . We define λ := maxn∈[N−M ] |λn|. We denote U by the eigenspace associated with λN and assume that U satisfies Assumption 1. Then, we have dM(X(l+1)) ≤ slλdM(X(l)). Remark 1. It is known that any Markov process on finite states converges to a unique distribution (equilibrium) if it is irreducible and aperiodic (see e.g., Norris (1998)). Theorem 1 includes this proposition as a special case with M = 1, C = 1, and Wl = 1 for all l ∈ N+. This is essentially the direct consequence of Perron – Frobenius’ theorem (see e.g., Meyer (2000)). See Appendix F.

heading: 4 APPLICATION TO GCN
text: We formulate a GCN (Kipf & Welling, 2017) without readout operations (Gilmer et al., 2017) using the dynamical system in the previous section and derive a sufficient condition in terms of the spectra of underlying graphs in which layer stacking nor non-linearity are not helpful for node classification.
Let G = (V,E) be an undirected graph where V is a set of nodes and E is a set of edges. We denote the number of nodes inG byN = |V | and identify V with [N ] by fixing an order of V . We associate a C dimensional signal to each node. X in the previous section corresponds to concatenation of the signals. GCNs iteratively update signals on V using the connection information and weights.
Let A := (1{(i,j)∈E})i,j∈[N ] ∈ RN×N be the adjacency matrix and D := diag(deg(i)i∈[N ]) ∈ RN×N be the degree matrix of G where deg(i) := |{j ∈ V | (i, j) ∈ E}| is the degree of node i. Let Ã := A + IN , D̃ := D + IN be the adjacent and degree matrix of graph G augmented with self-loops. We define the augmented normalized Laplacian (Wu et al., 2019a) of G by ∆̃ := IN−D̃− 1 2 ÃD̃− 1 2 and set P := IN−∆̃. Let L,C ∈ N+ be the layer and channel sizes, respectively. For weights Wl ∈ RC×C (l ∈ [L]), we define a GCN2 associated with G by f = fL ◦ · · · ◦ f1 where fl : RN×C → RN×C is defined by fl(X) := σ(PXWl). We are interested in the asymptotic behavior of the output X(L) of the GCN as L→∞. Suppose G has M connected components and let V = V1 t · · · t VM be the decomposition of the node set V into connected components. We denote an indicator vector of the m-th connected component by um := (1{n∈Vm})n∈[N ] ∈ RN . The following proposition shows that GCN satisfies the assumption of Corollay 3 (see Appendix B for proof). Proposition 1. Let λ1 ≤ · · · ≤ λN be the eigenvalue of P sorted in ascending order. Then, we have −1 < λ1, λN−M < 1, and λN−M+1 = · · · = λN = 1. In particular, we have λ := maxn=1,...,N−M |λn| < 1. Further, em := D̃ 1 2um for m ∈ [M ] are the basis of the eigenspace associated with the eigenvalue 1. Theorem 2. For any initial value X(0), the output of l-th layer X(l) satisfies dM(X(l)) ≤ (sλ)ldM(X (0)). In particular, dM(X(l)) exponentially converges to 0 when sλ < 1.
In the context of node classification tasks, we can interpret this corollary as the “information loss” of GCNs in the limit of infinite layers. For any X ∈M, if two nodes i, j ∈ V are in a same connected component and their degrees are identical, then, the column vectors of X that correspond to nodes i and j are identical. It means that we cannot distinguish these nodes using X . In this sense,M only has information about connected components and node degrees and we can interpret this theorem as the exponential information loss of GCNs in terms of the layer size. Similarly to the discussion in the previous section, X(l) converges to the trivial fixed point 0 when s < 1 (remember λN = 1). An interesting point is that even if s ≥ 1, X(l) can suffer from this information loss when s < λ−1. We note that the rate sλ in Theorem 2 depends on the spectra of the augmented normalized Laplacian, which is determined by the topology of the underlying graph G. Hence, our result explicitly relates the topological information of graphs and asymptotic behaviors of graph NNs. Remark 2. The old preprint (version 2) of Luan et al. (2019) formulated a theorem that explains the over-smoothing of non-linear GNNs. Specifically, it claimed that if a graph does not have a bipartite component and the input distribution is continuous, the rank of the output of a GCN converges to the number of connected components of the underlying graph as the layer size goes to infinity almost surely. However, it is not true in general as we give a counterexample in Appendix C.

heading: 5 ASYMPTOTIC BEHAVIOR OF GCN ON ERDŐS – RÉNYI GRAPH
text: Theorem 2 gives us a way to characterize the asymptotic behaviors of GCNs via the spectral distributions of the underlying graphs. To demonstrate this, we consider an Erdős – Rényi graph GN,p (Erdös & Rényi, 1959; Gilbert, 1959), which is a random graph that has N nodes and whose edges between two distinct nodes appear independently with probability p ∈ [0, 1], as an example. First, consider a (non-random) graph G with M connected components. Let 0 = µ̃1 = · · · = µ̃M < µ̃M+1 ≤ · · · ≤ µ̃N < 2 be eigenvalues of the augmented normalized Laplacian of G (see, Proposition 1) and set λ := minm=M+1,...,N |1 − µ̃m|(< 1). By Theorem 2, the output of GCN “loses information” as the layer size goes to infinity when the largest singular values of weights are strictly smaller than λ−1. Therefore, the closer the positive eigenvalues µm are to 1, the broader range of GCNs satisfies the assumption of Theorem 2.
2Following the original paper (Kipf & Welling, 2017), we use one-layer MLPs (i.e., Hl = 1 for all l ∈ N+.). However, our result holds for the multi-layer case
For an Erdős – Rényi graph GN,p, Chung & Radcliffe (2011) showed that when logNNp = o(1), the eigenvalues of the (usual) normalized Laplacian except for the smallest one converge to 1 with high probability (see Theorem 2 therein)3. We can interpret this theorem as the convergence of ErdősRényi graphs to the complete graph in terms of graph spectra. We can prove that the augmented normalized Laplacian behaves similarly (Lemma 6). By combining this fact with the discussion in the previous paragraph, we obtain the asymptotic behavior of GCNs on the Erdős – Rényi graph. See Appendix D for the complete proof.
Theorem 3. Consider a GCN on the Erdős-Rényi graph GN,p such that logNNp = o(1) as a function of N . For any ε > 0, if the supremum s of the maximum singular values of weights in the GCN
satisfies s < s0 := 17 √ Np−p+1 log(4N/ε) , then, for sufficiently large N , the GCN satisfies the condition of Theorem 2 with probability at least 1− ε.
Theorem 3 requires that an underlying graph is not extremely sparse. For example, suppose the node size is N = 20, 000, which is the approximately the maximum node size of datasets we use in experiments, and the edge probability is p = logN/N . Then, each node has the order of Np ≈ 4.3 adjacent nodes.
Under the condition of Theorem 3, the upper bound s0 →∞ as N →∞. It means that if the graph is sufficiently large and not extremely sparse, most GCNs suffer from the information loss. For the dependence on the edge probability p, s0 is an increasing function of p, which means the denser a graph is, the more quickly graph convolution operations mix signals on nodes and move them close to each other.
Theorem 3 implies that graph NNs perform poorly on dense NNs. More aggressively, we can hypothesize that the sparsity of practically available graphs is one of the reasons for the success of graph NNs in node classification tasks. To confirm this hypothesis, we artificially add edges to citation networks to make them dense in the experiments and observe the failure of graph NNs as expected (see Section 6.3).

heading: 6 EXPERIMENT

heading: 6.1 SYNTHESIS DATA: ONE-STEP TRANSITION
text: We numerically investigate how the transition f(X) := σ(PXW ) changes inputs using the vector field V (X) := f(X) − X4. For this purpose, we set N = 2, M = 1, and C = 1. Let λ1 ≤ λ2 be the eigenvalues of P . We choose W as |λ2|−1 ≤ W < |λ1|−1 so that Theorem 1 is applicable but is not reduced to the trivial situation (see, Appendix E.2). We choose the eigenvector e ∈ R2
3Chung et al. (2004) and Coja-Oghlan (2007) proved similar theorems. 4Since we consider the one-step transition only, we omit the subscript l from fl, Xl, and Wl.
associated with λ2 in two ways as described below. See Appendix H.1 for the concrete values of P , e, and W . Figure 1 shows the visualization of V . First, we choose the non-negative eigenvector e so that it satisfies Assumption 1 (Case 1). We see that the transition function f uniformly decreases the distance from M. This is consistent with the consequence of Theorem 1. Next, we choose the eigenvector e = [e1 e2]
> such that the signs of e1 and e2 differ (Case 2), which violates Assumption 1. We see that M is not invariant under f and f does not uniformly decrease the distance fromM. Therefore, we cannot remove the non-negativity assumption from Theorem 1.

heading: 6.2 SYNTHESIS DATA: DISTANCE TO INVARIANT SPACE
text: We evaluate the distance to the invariant spaceM using synthesis data. We randomly generate an Erdős – Rényi graph, a GCN on it, and an input signalX(0). We compute the distance between the lth intermediate output X(l) and the invariant spaceM for various edge probability p and maximum singular value s. Figure 2 plots the logarithm of the relative distance y(l) = log dM(X
(l)) dM(X(0))
with respect to the layer index l. From Theorem 1, we know that it is upper bounded by y(l) = l log(sλ). We see that this bound well approximates the actual value when sλ is small. On the other hand, it is loose for large sλ. We leave tighter bounds for dM in such a case for future research.

heading: 6.3 REAL DATA: EFFECT OF MAXIMUM SINGULAR VALUES ON PERFORMANCE
text: Theorem 2 implies that if s is smaller than the threshold λ−1, we cannot expect deep GCN to achieve good prediction accuracy. Conversely, if we can successfully train the model, s should avoid the region s ≤ λ−1. We empirically confirm these hypotheses using real datasets.
We use Cora, CiteSeer, and PubMed (Sen et al., 2008), which are standard citation network datasets. The task is to classify the genre of papers using word occurrences and citation relationships. We regard each paper as a node and citation relationship as an edge. Due to space constraints, we focus on Cora in the main article. See Appendix H.3 and I.3 for the other datasets. The discussion in Section 5 implies that Theorem 2 can support a wide range of GCNs when the underlying graph is relatively dense. However, the citation networks are too sparse to examine the aforementioned hypotheses — Theorem 2 gives a non-trivial result only when 1 ≤ s < λ−1 ≈ 1 + 3.62× 10−3. To circumvent this, we make noisy versions of citation networks by randomly adding edges to graphs. Through this manipulation, we can increase the value of λ−1 to 1.11.
Figure 3 (left) shows the accuracy for the test dataset in terms of the maximum singular values and the number of graph convolution layers. We can observe that when GCNs whose maximum singular value s is out of the region s < λ−1 outperform those inside the region in almost all configurations. Furthermore, the accuracy of GCNs with s = 10 are better than those without normalization (unnormalized). Figure 3 (right) shows the transition of the maximum singular values of the weights during training when we use a three-layered GCN. We can observe that the maximum singular value s does not shrink to the region s ≤ λ−1. In addition, when the layer size is small and predictive accuracy is high, GCNs gradually increase s from the initial value and avoid the region. In conclusion, the experiment results are consistent with the theorems.

heading: 6.4 REAL DATA: EFFECT OF SIGNAL COMPONENT PERPENDICULAR TO INVARIANT SPACE
text: We can decompose the output X of a model as X = X0 +X1 (X0 ∈ M, X1 ∈ M⊥). According to the theory, X0 has limited information for node classification. We hypothesize that the model emphasizes the perpendicular component X1 to perform good predictions. To quantitatively evaluate it, we define the relative magnitude of the perpendicular component of the output X by t(X) := X1/X0. Figure 4 compares this quantity and the prediction accuracy on the noisy version of Cora (see Appendix I.4 for other datasets). We observe that these two quantities are correlated (R = 0.545). If we remove GCNs have only one layer (corresponding to right points in the figure), the correlation coefficient is 0.827. This result does not contradict to the hypothesis above 5.

heading: 7 DISCUSSION
text: Applicability to Graph NNs on Sparse Graphs. We have theoretically and empirically shown that when the underlying graph is sufficiently dense and large, the threshold λ−1 is large (Theorem 2 and Section 6.3), which means many graph CNNs are eligible. However, real-world graphs are not often dense, which means that Theorem 2 is applicable to very limited GCNs. In addition, CojaOghlan (2007) theoretically proved that if the expected average degree of GN,p is bounded, the smallest positive eigenvalue of the normalized Laplacian of GN,p is o(1) with high probability. The asymptotic behaviors of graph NNs on sparse graphs are left for future research.
Remedy for Over-smoothing. Based on our theory, we can propose several techniques for mitigating the over-smoothing phenomena. One idea is to (randomly) sample edges in an underlying graph. The sparsity of practically available graphs could be a factor in the success of graph NNs. Assuming this hypothesis is correct, there is a possibility that we can relive the effect of over-smoothing by sparsification. Since we can never restore the information in pruned edges if we remove them permanently, random edge sampling could work better as FastGCN (Chen et al., 2018a) and GraphSAGE (Hamilton et al., 2017) do. Another idea is to scale node representations (i.e., intermediate or final output of graph NNs) appropriately so that they keep away from the invariant spaceM. Our proposed weight scaling mechanism takes this strategy. Recently, Zhao & Akoglu (2020) has pro-
5We cannot conclude that large perpendicular components are essential for good performance, since the maximum singular value s is correlated to the accuracy, too.
posed PairNorm to alleviate the over-smoothing phenomena. Although the scaling target is different – they rescaled signals whereas we normalized weights – theirs and ours have similar spirits.
Graph NNs with Large Weights. Our theory suggests that the maximum singular values of weights in a GCN should not be smaller than a threshold λ−1 because it suffers from information loss for node classification. On the other hand, if the scale of weights are very large, the model complexity of the function class represented by graph NNs increases, which may cause large generalization errors. Therefore, from a statistical learning theory perspective, we conjecture that the graph NNs with too-large weights perform poorly, too. A trade-off should exist between the expressive power and model complexity and there should be a “sweet spot” on the weight scale that balances the two.
Relation to Double Descent Phenomena. Belkin et al. (2019) pointed out that modern deep models often have double descent risk curves: when a model is under-parameterized, a classical biasvariance trade-off occurs. However, once the model has a large capacity and perfectly fits the training data, the test error decreases as we increase the number of parameters. To the best of our knowledge, no literature reported the double descent phenomena for graph NNs (it is consistent with the picture of the classical U-shaped risk curve in the previous paragraph). It is known that double descent phenomena do not occur in some situations, especially depending on regularization types. For example, while Belkin et al. (2019) employed the interpolating hypothesis with the minimum norm, Mei & Montanari (2019) found that the double descent was alleviated or disappeared when they used Ridge-type regularization techniques. Therefore, one can hypothesize the over-smoothing is a cause or consequence of regularization that is more like a Ridge-type rather than minimum-norm inductive bias.
Limitations in Graph NN Architectures. Our analysis is limited to graph NNs with the ReLU activation function because we implicitly use the property that ReLU is a projection onto the cone {X ≥ 0} (Appendix A, Lemma 3). This fact enables the ReLU function to get along with the nonnegativity of eigenvectors associated with the largest eigenvalues. Therefore, it is far from trivial to extend our results to other activation functions such as the sigmoid function or Leaky ReLU (Maas et al., 2013). Another point is that our formulation considers the update operation (Gilmer et al., 2017) of graph NNs only and does not take readout operations into account. In particular, we cannot directly apply our theory to graph classification tasks in which each sample is a graph.
Over-smoothing of Residual GNNs. Considering the correspondence of graph NNs and Markov processes (see Appendix F), one can imagine that residual links do not contribute to alleviating the over-smoothing phenomena because adding residual connections to a graph NN corresponds to converting a Markov process to its lazy version. When a Markov process converges to a stable distribution, the corresponding lazy process also converges eventually under certain conditions. It implies that residual links might not be helpful. However, Li et al. (2019) reported that graph NNs with as many as 56 layers performed well if they added residual connections. Considering that, the situation could be more complicated than our intuitions. The analysis of the role of residual connections in graph NNs is a promising direction for future research.

heading: 8 CONCLUSION
text: In this paper, to understand the empirically observed phenomena that deep non-linear graph NNs do not perform well, we analyzed their asymptotic behaviors by interpreting them as a dynamical system that includes GCN and Markov process as special cases. We gave theoretical conditions under which GCNs suffer from the information loss in the limit of infinite layers. Our theory directly related the expressive power of graph NNs and topological information of the underlying graphs via spectra of the Laplacian. It enabled us to leverage spectral and random graph theory to analyze the expressive power of graph NNs. To demonstrate this, we considered GCN on the Erdős – Rényi graph as an example and showed that when the underlying graph is sufficiently dense and large, a wide range of GCNs on the graph suffer from information loss. Based on the theory, we gave a principled guideline for how to determine the scale of weights of graph NNs and empirically showed that the weight normalization implied by our theory performed well in real datasets. One promising direction of research is to analyze the optimization and statistical properties such as the generalization power (Verma & Zhang, 2019) of graph NNs via spectral and random graph theories.

heading: ACKNOWLEDGMENTS
text: We thank Katsuhiko Ishiguro for providing a part of code for the experiments, Kohei Hayashi and Haru Negami Oono for giving us feedback and comments on the draft, Keyulu Xu and anonymous reviewers for fruitful discussions via OpenReview, and Ryuta Osawa for pointing out errors and suggesting improvements of the paper. TS was partially supported by JSPS KAKENHI (15H05707, 18K19793, and 18H03201), Japan Digital Design, and JST CREST.

heading: A PROOF OF THEOREM 1
text: As we wrote in the main article, it is enough to show the following lemmas (definition of miscellaneous variables are as in Section 3.2). Remember that λ = supn∈[N−M ] |λn| and slh is the maximum singular value of Wlh Lemma 1. For any X ∈ RN×C , we have dM(PX) ≤ λdM(X). Lemma 2. For any X ∈ RN×C , we have dM(XWlh) ≤ slhdM(X). Lemma 3. For any X ∈ RN×C , we have dM(σ(X)) ≤ dM(X).
Proof of Lemma 1. Since P is a symmetric linear operator on U⊥, we can choose the orthonormal basis (em)m=M+1,...,N of U⊥ consisting of the eigenvalue of P |U⊥ . Let λm be the eigenvalue of P to which em is associated (m = M + 1, . . . , N ). Note that since the operator norm of P |U⊥ is λ, we have |λm| ≤ λ for all m = M + 1, . . . , N . Since (em)m∈[N ] forms the orthonormal basis of
RN , we can uniquely write X ∈ RN×C as X = ∑N m=1 em ⊗ wm for some wm ∈ RC . Then, we
have d2M(X) = ∑N m=M+1 ‖wm‖2 where ‖ · ‖ is the 2-norm of a vector. On the other hand, we have
PX = N∑ m=1 Pem ⊗ wm
= M∑ m=1 Pem ⊗ wm + N∑ m=M+1 Pem ⊗ wm
= M∑ m=1 Pem ⊗ wm + N∑ m=M+1 em ⊗ (λmwm)
Since U is invariant under P , for any m ∈ [M ], we can write Pem as a linear combination of en(n ∈ [M ]). Therefore, we have d2M(PX) = ∑N m=M+1 ‖λmwm‖2. Then, we obtain the desired inequality as follows:
d2M(PX) = N∑ m=M+1 ‖λmwm‖2
≤ λ2 N∑
m=M+1
‖wm‖2
≤ λ2 N∑
m=M+1
‖wm‖2
= λ2d2M(X).
Proof of Lemma 2. Using the same decomposition of X as the proof in Lemma 1, we have
XWlh = N∑ m=1 em ⊗ (W>lhwm)
= M∑ m=1 em ⊗ (W>lhwm) + N∑ m=M+1 em ⊗ (W>lhwm).
Therefore, we have
d2M(XWlh) = N∑ m=M+1 ‖W>lhwm‖2
≤ s2lh N∑
m=M+1
‖wm‖2
= s2lhd 2 M(X).
Proof of Lemma 3. We choose (em)m=N−M+1,...,N as in the proof of Lemma 1. We denote X = (Xnc)n∈[N ],c∈[C] and en = (emn)m∈[N ], respectively. Let (e′c)c∈[C] be the standard basis of RC . Then, (en ⊗ e′c)n∈[N ],c∈[C] is the orthonormal basis of RN×C , endowed with the standard inner product as a Euclid space. Therefore, we can decompose X as X = ∑N n=1 ∑C c=1 ancen⊗ e′c where
anc = 〈X, en ⊗ e′c〉 = ∑N m=1Xmcemn. Then, we have d 2 M(X) = ∑N n=M+1 ‖ ∑C c=1 ance ′ c‖2,
which is further transformed as
d2M(X) = N∑ n=M+1 ∥∥∥∥∥ C∑ c=1 ance ′ c ∥∥∥∥∥ 2
= N∑ n=M+1 C∑ c=1 a2nc
= C∑ c=1 ( N∑ n=1 a2nc − M∑ n=1 a2nc )
= C∑ c=1
( ‖X·c‖2 −
M∑ n=1
〈X·c, en〉2 ) ,
where X·c is the c-th column vector of X . Similarly, we have
d2M(σ(X)) = C∑ c=1 ( ‖X+·c‖2 − M∑ n=1 〈X+·c , en〉2 ) ,
where we denote σ(X) = (X+nc)n∈[N ],c∈[C] in shorthand. Therefore, the inequality follow from the following lemma.
Lemma 4. Let x ∈ RN and v1, . . . , vM ∈ RN be orthonormal vectors (i.e., 〈vm, vn〉 = δmn) satisfying vm ≥ 0 for allm ∈ [M ]. Then, we have ‖x‖2− ∑M m=1〈x, vm〉2 ≥ ‖x+‖2− ∑M m=1〈x+, vm〉2 where x+ := max(x, 0) for x ∈ R. Proof. The value ‖y‖2− ∑M m=1〈y, um〉2 is invariant under simultaneous coordinate permutation of y and um’s. Therefore, we can assume without loss of generality that the coordinate of x are sorted: x1 ≤ . . . ≤ xL < 0 ≤ xL+1 ≤ · · · ≤ xN for some L ≤ N . Then, we have
‖x‖2 − ‖x+‖2 = L∑ n=1 x2n. (1)
When L = 0, the sum in the right hand side is treated as 0. On the other hand, writing as vm = (vnm)n∈[N ], direct calculation shows
M∑ m=1 〈x, vm〉2 − 〈x+, vm〉2 = M∑ m=1 ( L∑ n=1 xnvnm )2 − 2 L∑ n=1 N∑ l=L+1 xnxlvnmvlm  . (2) Let Im := {n ∈ [N ] | vnm > 0} be the support of vm form ∈ [M ]. We note that ifm 6= m′ ∈ [M ], we have Im ∩ Im′ = ∅ since if there existed n ∈ Im ∩ Im′ , we have
0 = 〈vm, vm′〉 ≥ vnmvnm′ > 0, which is contradictory. Therefore, N∑ m=1 ( L∑ n=1 xnvnm )2 = N∑ m=1  ∑ n∈Im∩[L] xnvnm 2
≤ N∑ m=1  ∑ n∈Im∩[L] x2n  ∑ n∈Im∩[L] v2nm  (∵ Cauchy–Schwarz inequality) ≤
N∑ m=1  ∑ n∈Im∩[L] x2n  (∵ ‖vm‖2 = 1) ≤
L∑ n=1 x2n. (3)
We used the fact that Im’s are disjoint and vnm = 0 if n 6∈ ∪mIm in the first equality above. Further, we have xnxlvnmvlm ≤ 0 for 1 ≤ n ≤ L and L + 1 ≤ l ≤ N by the definition of L and non-negativity of vm. By combining (1), (2), and (3), we have
M∑ m=1 〈x, vm〉2 − 〈x+, vm〉2 ≤ L∑ n=1 x2n = ‖x‖2 − ‖x+‖2.
Proof of Theorem 1. By Lemma 1, 2, and 3, we have
dM(fl(X)) = dM(σ(· · ·σ(σ(︸ ︷︷ ︸ H times PX)Wl1)Wl2 · · ·WlHl))
≤ dM(σ(· · ·σ(σ(︸ ︷︷ ︸ H−1 times PX)Wl1)Wl2 · · · )WlHl) ≤ slHl−1dM(σ(· · ·σ(σ(︸ ︷︷ ︸ H−1 times PX)Wl1)Wl2 · · · )WlHl−1))
· · ·
≤ ( Hl∏ h=1 slh ) dM(PX)
≤ sldM(PX) ≤ slλdM(X).

heading: B PROOF OF PROPOSITION 1
text: Proof. Let µ̃1 ≤ · · · ≤ µ̃N be the eigenvalue of the augmented normalized Laplacian ∆̃, sorted in ascending order. Since P = IN − ∆̃, it is enough to show µ̃1 = · · · = µ̃M = 0, µ̃M+1 > 0, and µ̃N < 2. For the first two, the statements are equivalent to that ∆̃ is positive semi-definite and that the multiplicity of the eigenvalue 0 is same as the number of connected components 6. This is well-known for Laplacian or its normalized version (see, e.g., Chung & Graham (1997)) and the proof for ∆̃ is similar. By direct calculation, we have
x>∆̃x = 1
2 N∑ i,j=1 aij ( xi√ di + 1 − xj√ dj + 1 )2
for any x = [x1 · · · xN ]> ∈ RN . Therefore, ∆̃ is positive semi-definite and hence µ̃1 ≥ 0.
Suppose temporally that G is connected. If x ∈ RN is an eigenvector associated to 0, then, by the aformentioned calculation, xi√
di+1 and xj√ dj+1 must be same for all pairs (i, j) such that aij >
0. However, since G is connected, xi√ di+1 must be same value for all i ∈ [N ]. That means the multiplicity of the eigenvalue 0 is 1 and any eigenvector associated to 0 must be proportional to D̃ 1 21. Now, suppose G has M connected components V1, . . . , VM . Let ∆̃m be the augmented normalized Laplacians corresponding to each connected component Vm for m ∈ [M ]. By the aformentioned discussion, ∆̃m has the eigenvalue 0 with multiplicity 1. Since ∆̃ is the direct sum of ∆̃′ms, the eigenvalue of ∆̃ is the union of those for ∆̃m’s. Therefore, ∆̃ has the eigenvalue 0 with multiplicity M and em = D̃ 1 21m’s are the orthogonal basis of the eigenspace.
6The former statement is identical to Lemma 1 and latter one is the extension of Lemma 2 of Wu et al. (2019a).
Finally, we prove µ̃N < 2. Let µN be the largest eigenvalue of the normalized Laplacian ∆ = D− 1 2 (D −A)D− 12 , where D− 12 ∈ RN×N is the diagonal matrix defined by
D − 12 ii =
{ deg(i)− 1 2 (if deg(i) 6= 0)
0 (if deg(i) = 0) .
Note that D− 1 2D 1 2 nor D 1 2D− 1 2 are not equal to the identity matrix IN in general. However, we have
L = D 1 2D− 1 2LD− 1 2D 1 2 (4)
where L = D −A is the (unnormalized) Laplacian. Therefore, we have
µ̃N = max x 6=0
x>∆̃x
‖x‖
= max x6=0
x>D̃− 1 2LD̃− 1 2x
‖x‖
= max x6=0
x>D̃− 1 2D 1 2D− 1 2LD− 1 2D 1 2 D̃− 1 2x
‖x‖ (∵ (4))
= max x6=0
(D 1 2 D̃− 1 2x)>∆(D 1 2 D̃− 1 2x)
‖x‖
= max x 6=0
D 1 2 D̃− 1 2 x6=0
(D 1 2 D̃− 1 2x)>∆(D 1 2 D̃− 1 2x)
‖x‖
= max x 6=0
D 1 2 D̃− 1 2 x6=0
(D 1 2 D̃− 1 2x)>∆(D 1 2 D̃− 1 2x) ‖D 12 D̃− 12x‖ ‖D 12 D̃− 12x‖ ‖x‖
≤ max x 6=0
D 1 2 D̃− 1 2 x6=0
(D 1 2 D̃− 1 2x)>∆(D 1 2 D̃− 1 2x)
‖D 12 D̃− 12x‖ max x 6=0
D 1 2 D̃− 1 2 x 6=0
‖D 12 D̃− 12x‖ ‖x‖
≤ max y 6=0
y>∆y
‖y‖ max x 6=0 ‖D 12 D̃− 12x‖ ‖x‖
= µN max n∈[N ]
( deg(i)
deg(i) + 1
) 1 2
≤ µN .
Therefore, we have µ̃N ≤ µN 7. Since maxi∈[N ] ( deg(i) deg(i)+1 ) 1 2
< 1, the equality µ̃N = µN holds if and only if µN = 0, that is, G has N connected components. On the other hand, it is known that µN ≤ 2 and the equality holds if and only if G has non-trivial bipartite graph as a connected component (see, e.g., Chung & Graham (1997)). Therefore, µ̃N = µN and µN = 2 does not hold simultaneously and we obtain µN < 2.

heading: C COUNTEREXAMPLE OF PREVIOUS STUDY ON OVER-SMOOTHING FOR NON-LINEAR GNNS
text: We restate Theorem 1 of the preprint (version2) of Luan et al. (2019)8. LetG be a simple undirected graph with N nodes and k connected components such that it does not have a bipartite component. Let L = D̃−1/2ÃD̃−1/2 ∈ RN×N be the augmented normalized Laplacian of G. Let F ∈ N+ and
7Theorem 1 of Wu et al. (2019a) showed that this inequality strictly holds when G is simple and connected. We do not require this assumption.
8https://arxiv.org/abs/1906.02174v2
Wn ∈ RF×F be the weight of the n-th layer for n ∈ N+. For the input X ∈ RN×F , we define the output Yn ∈ RN×F of the n-th layer of a GCN by Yn = σ(L · · ·σ(LXW0) · · ·Wn) where σ is the ReLU function. We assume the input X is drawn from a continuous distribution on RN×F . Then, the theorem claims that we have limn→∞ rank(Yn) = k almost surely with respect to the distribution of X .
We construct a conterexample. Consider a graph G consisting of N = 4 nodes whose adjacency matrix is
A = 1 1 1 11 1 1 01 1 1 0 1 0 0 1  . Note that G is connected (i.e., k = 1) and is not bipartite. We make a GCN with F = 3 channels and whose weight matrices are Wn = I3 (the identity matrix of size 3) for all n ∈ N. For the distribution of the input X , we consider an absolutely continuous distribution with respect to the Lebesgue measure on R4×3 such that P (X ≥ 0) > 0 (here, X ≥ 0 means the element-wise comparison). For example, the standard Gaussian distribution satisfies the condition.
Since L ≥ 0, we have Yn = LnX if X ≥ 0. Let L = P>ΛP be the diagonalization of L where P ∈ O(4) is an orthogonal matrix of size 4. Since rank(L) = 3, we have rank(Λn) = 3 for any n (we can assume that Λ44 = 0 without loss of generality). Therefore, under the condition X ≥ 0, we have
rank(Yn) = 3 ⇐⇒ rank(P>ΛnPX) = 3
⇐⇒ X ∈ {P−1 [B v]> | B ∈ R3×3 is invertible, v ∈ R3}.
Note that the last condition is independent of n. Since the set of invertible matrices is dense in the set of all matrices of the same size (with respect to the standard topology of the Euclidean space), we have P ({rank(Yn) = 3 for all n ∈ N}) > 0. Therefore, we have limn→∞ rank(Yn) = 3 with a non-zero probability.

heading: D PROOF OF THEOREM 3
text: We follow the proof of Theorem 2 of Chung & Radcliffe (2011). The idea is to relate the spectral distribution of the normalized Laplacian with that of its expected version. Since we can compute the latter one explicitly for the Erdős-Rényi graph, we can derive the convergence of spectra. We employ this technique and derive similar conclusion for the augmented normalized Laplacian.
First, we consider genral random graphs not restricted to Erdős-Rényi graphs. Let N ∈ N+, and P = (pij)i,j∈[N ] be a non-negative symmetric matrix (meaning that pij ≥ 0 for any i, j ∈ [N ]). Let G be an undirected random graph with N nodes such that an edge between i and j is independently present with probability pij . Let A and D be the adjacency and the degree matrices of G, respectively (that is, Aij ∼ Ber(pij), i.i.d.). Define the expected node degree of node i by ti := ∑N j=1 pij . Let Ã := A + IN , D̃ := D + IN and define Ā := E[Ã] = P + IN and D̄ := E[D̃] = diag(t1, . . . , tN ) + IN correspondingly. We define the augmented normalized Laplacian ∆̃ of G by ∆̃ := IN − D̃− 1 2 ÃD̃− 1 2 and its expected version by ∆̄ := IN − D̄− 1 2 ĀD̄− 1 2 9. For a symmetric matrix X ∈ RN , we define its eigenvalues, sorted in ascending order by λ1(X) ≤ · · · ≤ λN (X) and its operator norm by ‖X‖ = maxn∈[N ] |λn(X)|. Lemma 5 (Ref. Chung & Radcliffe (2011) Theorem 2). Let δ := minn∈[N ] tn be the minimum expected degree of G. Set k(ε) := 3(1 + log(4/ε)). Then, for any ε > 0, if δ + 1 > k(ε) logN , we have
max n∈[N ] ∣∣∣λn(∆̃)− λn(∆̄)∣∣∣ ≤ 4√3 log(4N/ε) δ + 1
with probability at least 1− ε. 9Note that E[∆̃] 6= ∆̄ in general due to the dependence between Ã and D̃.
Proof. By Weyl’s theorem, we have maxn∈[N ] ∣∣∣λn(∆̃)− λn(∆̄)∣∣∣ ≤ ‖∆̃ − ∆̄‖. Therefore, it is
enough to bound ‖∆̃ − ∆̄‖. Let C := IN − D̄− 1 2 ÃD̄. By the triangular inequality, we have ‖∆̃− ∆̄‖ ≤ ‖∆̃− C‖+ ‖C − ∆̄‖. We will bound these terms respectively.
First, we bound ‖C−∆̄‖. Direct calculation showsC−∆̄ = −D̄− 12 (A−P )D̄− 12 . LetEij ∈ RN×N be a matrix defined by
(Eij)kl = { 1 if (i = k and i = l) or (i = l and j = k), 0 otherwise.
We define the random variable Yij by
Yij := Aij − pij√ ti + 1 √ tj + 1 Eij .
Then, Yij’s are independent and we have C − ∆̄ = ∑N i,j=1 Yij . To apply Theorem 5 of Chung &
Radcliffe (2011) to Yij’s, we bound ‖Yij − E[Yij ]‖ and ‖ ∑N i,j=1 E[Y 2ij ]‖. First, we have
‖Yij − E[Yij ]‖ = ‖Yij‖ ≤ ‖Eij‖√ ti + 1 √ tj + 1 ≤ (δ + 1)−1.
Since
E[Y 2ij ] = pij − p2ij
(ti + 1)(tj + 1)
{ Eii + Ejj (if i 6= j), Eii (if i = j),
we have ∥∥∥∥∥∥ N∑
i,j=1
E[Y 2ij ] ∥∥∥∥∥∥ = ∥∥∥∥∥∥ N∑ i,j=1 pij − p2ij (ti + 1)(tj + 1) Eii ∥∥∥∥∥∥ = max i∈[N ]  N∑ j=1 pij − p2ij (ti + 1)(tj + 1)
 ≤ max i∈[N ]  N∑ j=1 pij (ti + 1)(tj + 1)
 ≤ (δ + 1)−1.
By letting a← √
3 log(4N/ε) δ+1 , M ← (δ + 1) −1, v2 ← (δ + 1)−1 and applying Theorem 5 of Chung & Radcliffe (2011), we have
Pr(‖C − ∆̄‖ > a) ≤ 2N exp ( − a 2
2(δ + 1)−1 + 2(δ + 1)−1a/3 ) ≤ 2N exp ( −3 log(4N/ε)
2(1 + a/3)
) .
By the definition of k(ε), we have a < 1 if δ + 1 > k(ε) log n. For such δ, we have Pr(‖C − ∆̄‖ > a) ≤ 2N exp ( −3 log(4N/ε)
2(1 + a/3) ) ≤ 2N exp (− log(4N/ε)) (∵ a < 1)
= ε
2 . (5)
Next, we bound ‖∆̃ − C‖. First, since a < 1, by Chernoff bound (see, e.g. Angluin & Valiant (1979); Hagerup & Rüb (1990))), we have
Pr(|di − ti| > a(ti + 1)) ≤ 2 exp ( −a 2(ti + 1)
3 ) ≤ 2 exp ( −a 2(δ + 1)
3 ) = ε
2N .
Therefore, if |di − ti| ≤ a(ti + 1), then we have∣∣∣∣∣ √ di + 1 ti + 1 − 1 ∣∣∣∣∣ ≤ ∣∣∣∣di + 1ti + 1 − 1 ∣∣∣∣ (∵ |√x− 1| ≤ |x− 1| for x ≥ 0) =
∣∣∣∣di − titi + 1 ∣∣∣∣
≤ a. Therefore, by union bound, we have
‖D̄− 12 D̃ 12 − IN‖ = max i∈[N ] ∣∣∣∣∣ √ di + 1 ti + 1 − 1 ∣∣∣∣∣ ≤ a with probability at least 1− ε/2. Further, since the eigenvalue of the augmented normalized Laplacian is in [0, 2] by the proof of Proposition 1, we have ‖IN − ∆̃‖ ≤ 1. By combining them, we have
‖∆̃− C‖ = ‖(D̄− 12 D̃ 12 − IN )(IN − ∆̃)D̃ 1 2 D̄− 1 2 + (IN − ∆̃)(I − D̃ 1 2 D̄− 1 2 )‖
≤ ‖(D̄− 12 D̃ 12 − IN‖‖D̃ 1 2 D̄− 1 2 ‖+ ‖I − D̃ 12 D̄− 12 ‖
≤ a(a+ 1) + a. (6) From (5) and (6), we have
‖∆̃− ∆̄‖ ≤ ‖∆̃− C‖+ ‖C − ∆̄‖ ≤ a+ a(a+ 1) + a ≤ a2 + 3a ≤ 4a (∵ a < 1)
with probability at least 1− ε by union bound.
LetN ∈ N+ and p > 0. In the case of the Erdős-Rényi graphGN,p, we should set P = p(JN −IN ) where JN ∈ RN×N are the all-one matrix. Then, we have Ā = pJN + (1 − p)IN , D̄ = (Np − p + 1)IN , and ∆̄ = pNp−p+1 (NIN − JN ). Since the eigenvalue of JN is N (with multiplicity 1) and 0 (with multiplicity N − 1), the eigenvalue of ∆̄ is 0 (with multiplicity 1) and NpNp−p+1 (with multiplicity N − 1). For GN,p, δ is the expected average degree (N − 1)p. Hence, we have the following lemma from Lemma 5:
Lemma 6. Let ∆̃ be its augmented normalized Laplacian of the Erdős-Rényi graph GN,p. For any ε > 0, if Np−p+1logN > k(ε) := 3(1 + log(4/ε)), then, with probability at least 1− ε, we have
max i=2,...,N ∣∣∣∣λi(∆̃)− NpNp− p+ 1 ∣∣∣∣ ≤ 4 √ 3 log(4N/ε) Np− p+ 1 .
Corollary 4. Consider GCN on GN,p. Let Wl be the weight of the l-th layer of GCN and sl be the maximum singular value of Wl for l ∈ N+. Set s := supl∈N+ . Let ε > 0. We define k(ε) :=
3(1 + log(4/ε)) and l(N, p, ε) = 1−pNp−p+1 + 4 √ 3 log(4N/ε) Np−p+1 . If Np−p+1 logN > k(ε) and s ≤ l(N, ε)
−1, then, GCN on GN,p satisfies the assumption of Theorem 2 with probability at least 1− ε.
Proof of Theorem 3. Since logNNp = o(1), for fixed ε, we have
Np− p+ 1 logN > Np logN > k(ε)
for sufficiently large N . Further, Np→∞ as N →∞ when logNNp = o(1). Therefore, we have
(1− p)2
Np− p+ 1 ≤ 1 Np ≤ (7− 4
√ 3)2 log
( 4N
ε ) for sufficiently large N . Hence.
1− p Np− p+ 1 ≤ (7− 4 √ 3)
√ log(4N/ε)
Np− p+ 1 .
Therefore, we have l(N, p, ε) ≤ 7 √
log(4N/ε) Np−p+1 . Therefore, if s ≤ 1 7 √ Np−p+1 log(4N/ε) , then we have
s ≤ l(N, p, ε)−1.

heading: E MISCELLANEOUS PROPOSITIONS
text: E.1 INVARIANCE OF ORTHOGONAL COMPLEMENT SPACE
Proposition 2. Let P ∈ RN×N be a symmetric matrix, treated as a linear operator P : RN → RN . If a subspace U ⊂ RN is invariant under P (i.e., if u ∈ U , then Pu ∈ U ), then, U⊥ is invariant under P , too.
Proof. For any u ∈ U⊥ and v ∈ U , by symmetry of P , we have 〈Pu, v〉 = (Pu)>v = u>P>v = u>Pv = 〈u, Pv〉.
Since U is an invariant space of P , we have Pv ∈ U . Hence, we have 〈u, Pv〉 = 0 because u ∈ U⊥. We obtain Pu ∈ U⊥ by the definition of U⊥.

heading: E.2 CONVERGENCE TO TRIVIAL FIXED POINT
text: Let P ∈ RN×N be a symmetric matrix, Wl ∈ RC×C , sl be the maximum singular value of Wl for l ∈ N+. We define fl : RN×C → RN×C by fl(X) := σ(PXWl) where σ is the element-wise ReLU function. Proposition 3. Suppose further that the operator norm of P is no larger than λ, then we have ‖fl(X)‖F ≤ slλ‖X‖F for any l ∈ N+. In particular, let s := supl∈N+ sl. If sλ < 1, then, Xl exponentially approaches 0 as l→∞.
Proof. Since λ is the operator norm of P |U⊥ , the assumption implies that the operator norm of P itself is no larger than λ. Therefore, we have ‖PXWl‖F ≤ λ‖XWl‖F ≤ slλ‖X‖F. On the other hand, since σ(x)2 ≤ x2 for any x ∈ R, we have ‖σ(X)‖F ≤ ‖X‖F for any X ∈ RN×C . Combining the two, we have ‖fl(X)‖F ≤ ‖PXWl‖F ≤ slλ‖X‖F.

heading: E.3 STRICTNESS OF THEOREM 1
text: Theorem 1 implies that if sλ ≤ 1, then, one-step transition fl does not increase the distance toM. In this section, we first prove that this theorem is strict in the sense that, there exists a situation in which slλ > 1 holds and the distance dM increases by one-step transition fl at some point X .
Set N ← 2, C ← 1, and M ← 1 in Section 3.2. For µ, λ > 0, we set P ← [ µ 0 0 λ ] , e← [ 1 0 ] , U ← {[ x y ] | y = 0 } .
Then, by definition, we can check that the 3-tuple (P, e, U) satisfies the Assumptions 1 and 2. Set M := U ⊗ R = U and choose W ∈ R so that W > λ−1. Finally define f : RN×C → RN×C by f(X) := σ(PXW ) where σ is the element-wise ReLU function.
Proposition 4. We have dM(f(X)) > dM(X) for any X = [x1 x2] > ∈ R2 such that x2 > 0.
Proof. By definition, we have dM(X) = |x2|. On the other hand, direct calculation shows that fl(X) = [ (WµX1) + (WλX2) + ]>
and dM(fl(X)) = (WλX2)+ where x+ := max(x, 0) for x ∈ R. Since W > λ−1 and x2 > 0, we have dM(fl(X)) > dM(X).
Next, we prove the non-strictness of Theorem 1 in the sense that there exists a situation in which slλ > 1 holds and the distance dM uniformly decreases by fl. Again, we set Set N ← 2, C ← 1, and M ← 1. Let λ ∈ (1, 2) and set
P ← λ 2 [ 1 −1 −1 1 ] , e← 1√
2
[ 1 1 ] , U ← {[ x y ] | x = y } Then, we can directly show that 3-tuple (P, e, U) satisfies the Assumptions 1 and 2. Set W ← 1. Proposition 5. We have Wλ > 1 and dM(fl(X)) < dM(X) for all X ∈ R2.
Proof. First, note that e′ := 1√ 2 [1 −1]> is the eigenvector of P associated to λ: Pe′ = λe′. For X = ae + be′ (a, b ¿ 0), the distance between X andM is dM(X) = |b|. On the other hand, by direct computation, we have
f(X) = σ(PXW ) =  [ 0 λb√ 2 ]> (if b ≥ 0),[
−λb√ 2
0 ]> (if b < 0).
Therefore, the distance between f(X) and M is dM(f(X)) = λ|b|/2. Since λ < 2, we have dM(f(X)) < dM(X) for any X ∈ R2.
We have shown that the non-negativity of e (Assumption 1) is not a redundant condition in Section 6.1.

heading: F RELATION TO MARKOV PROCESS
text: It is known that any Markov process on finite states converges to a unique distribution (equilibrium) if it is irreducible and aperiodic (see, e.g., Norris (1998)). As we see in this section, this theorem is the special case of Corollary 3.
Let S := {1, . . . , N} be a finite discrete state space. Consider a Markov process on S characterized by a symmetric transition matrix P = (pij)i,j∈[N ] ∈ RN×N such that P ≥ 0 and P1 = 1 where 1 is the all-one vector. We interpret pij as the transition probability from a state i to j. We associate P with a graph GP = (VP , EP ) by VP = [N ] and (i, j) ∈ EP if and only if pij > 0. Since P is symmetric, we can regard GP as an undirected graph. We assume P is irreducible and aperiodic 10. Perron – Frobenius’ theorem (see, e.g., Meyer (2000)) implies that P satisfy the assumption of Corollary 3 with M = 1. Proposition 6 (Perron – Frobenius). Let the eigenvalues of P be λ1 ≤ · · · ≤ λN . Then, we have −1 < λ1, λN−1 < 1, and λN = 1. Further, there exists unique vector e ∈ RN such that e ≥ 0, ‖e‖ = 1, and e is the eigenvector for the eivenvalue 1. Corollary 5. Let λ := maxn=1,...,N−1 |λn|(< 1) andM := {e ⊗ w | w ∈ RC}. If sλ < 1, then, for any initial valueX1, Xl exponentially approachesM as l→∞.
If we set C = 1 and Wl = 1 for all l ∈ N+, then, we can inductively show that Xl ≥ 0 for any l ≥ 2. Therefore, we can interpret Xl as a measure on S. Suppose further that we take the initial value X1 as X1 ≥ 0 and X>1 1 = 1 so that we can interpret X1 as a probability distribution on S. Then, we can inductively show that Xl ≥ 0, X>l 1 = 1 (i.e., Xl is a probability distribution on
10A symmetric matrix A is called irreducible if and only if GA is connected. We say a graph G is aperiodic if the greatest common divisor of length of all loops in G is 1. A symmetric matrix A is aperiodic if the graph GA induced by A is aperiodic.
S), and Xl+1 = σ(PXlWl) = PXl for all l ∈ N+. In conclusion, the corollary is reduced to the fact that if a finite and discrete Markov process is irreducible and aperiodic, any initial probability distribution converges exponentially to an equibrilium. In addition, the the rate λ corresponds to the mixing time of the Markov process.

heading: G GCN DEFINED BY NORMALIZED LAPLACIAN
text: In Section 4, we defined P using the augmented normalized Laplacian ∆̃ by P = IN − ∆̃. We can alternatively use the usual normalized Laplacian ∆ instead of the augmented one to define P and want to apply the theory developed in Section 3.2. We write the normalized Laplacian version as P∆ := IN − ∆. The only obstacle is that the smallest eigenvalue λ1 of P∆ can be equal to −1, while that of P is strictly larger than −1 (see, Proposition 1). This corresponds to that fact the largest eigenvalue of ∆̃ is strictly smaller than 2, while that for ∆ can be 2. It is known that the largest eigenvalue of ∆ is 2 if and only if the graph has a non-trivial bipartite connected component (see, e.g., Chung & Graham (1997)). Therefore, we can develop a theory using the normalized Laplacian instead of the augmented one in parallel for such a graph G.
In Section 5, we characterized the asymptotic behavior of GCN defined by the augmented normalized Laplacian via its spectral distribution (Lemma 6 of Appendix D). We can derive a similar claim for GCN defined via the normalized Laplacian using the original theorem for the normalized Laplacian in Chung & Radcliffe (2011) (Theorem 7 therein). The normalized Laplacian version of GCN is advantegeous over the one made from the augmented one because we know its spectral distribution for broader range of random graphs. For example, Chung & Radcliffe (2011) proved the convergence of the spectral distribution of the normalized Laplacian for Chung-Lu’s model (Chung & Lu, 2002), which includes power law graphs as a special case (see, Theorem 4 of Chung & Radcliffe (2011)).

heading: H DETAILS OF EXPERIMENT SETTINGS

heading: H.1 EXPERIMENT OF SECTION 6.1
text: We set the eigenvalue of P to λ1 = 0.5 and λ2 = 1.0 and randomly generated P until the eigenvector e associated to λ2 satisfies the condition of each case described in the main article. We set W = 1.2 and used the following values for each case as P and e.
H.1.1 CASE 1
P = [ 0.7469915 0.2499819 0.2499819 0.7530085 ] , e = [ 0.7028392 −0.71134876 ] .
H.1.2 CASE 2
P = [ 0.6899574 −0.2426827 −0.2426827 0.8100426 ] , e = [ 0.61637234 −0.78745485 ] .

heading: H.2 EXPERIMENT OF SECTION 6.2
text: We randomly generated an Erdős – Rényi graph GN,p with N = 1000 and randomly generated a one-of-K hot vector for each node and embed it to a C-dimensional vector using a random matrix whose elements were randomly sampled from the standard Gaussian distribution. Here, K = 10 and C = 32. We treated the resulting single as the input signal X(0) ∈ RN×C . We constructed a GCN with L = 10 layers and C channels. All parameters were i.i.d. sampled from the Gaussian distribution whose standard deviation is same as the one used in LeCun et al. (2012)11 and multiplied a scalar to each weight matrix so that the largest singular value equals to a specified value s. We used three configurations (p, s) = (0.1, 0.1), (0.5, 1.0), (0.5, 10.0). λ of the generated GCNs are 0.063, 0.197, 0.194, respectively. See Appendix 6.2 for the results of other configurations of (p, s).
11This is the default initialization method for weight matrices in Chainer and Chainer Chemistry.

heading: H.3 EXPERIMENT OF SECTION 6.3

heading: H.3.1 DATASET
text: We used the Cora (McCallum et al., 2000; Sen et al., 2008), CiteSeer (Giles et al., 1998; Sen et al., 2008), and PubMed(Sen et al., 2008) datasets for experiments. We obtained the preprocessed dataset from the code repository of Kipf & Welling (2017)12. Table 1 summarizes specifications of datasets and their noisy version (explained in the next section).
The Cora dataset is a citation network dataset consisting of 2708 papers and 5429 links. Each paper is represented as the occurence of 1433 unique words and is associated to one of 7 genres (Case Based, Genetic Algorithms, Neural Networks, Probabilistic Methods, Reinforcement Learning, Rule Learning, Theory). The graph made from the citation links has 78 connected components and the smallest positive eigenvalue of the augmented Normalized Laplacian is approximately µ̃ = 3.62 × 10−3. Therefore, the upper bound of Theorem 2 is λ−1 = (1 − µ̃)−1 ≈ 1 + 3.62 × 10−3. 818 out of 2708 samples are labelled as “Probabilistic Methods”, which is the largest proportion. Therefore, the chance rate is 818/2708 = 30.2%.
The CiteSeer dataset is a citation network dataset consisting of 3312 papers and 4732 links. Each paper is represented as the occurence of 3703 unique words and is associated to one of 6 genres (Agents, AI, DB, IR, ML, HCI). The graph made from the citation links has 438 connected components and the smallest positive eigenvalue of the augmented Normalized Laplacian is approximately µ̃ = 1.25×10−3. Therefore, the upper bound of Theorem 2 is λ−1 = (1− µ̃)−1 ≈ 1+1.25×10−3. 701 out of 2708 samples are labelled as “IR”, which is the largest proportion. Therefore, the chance rate is 701/3312 = 21.1%.
The PubMed dataset is a citation network dataset consisting of 19717 papers and 44338 links. Each paper is represented as the occurence of 500 unique words and is associated to one of 3 genres (“Diabetes Mellitus, Experimental”, “Diabetes Mellitus Type 1”, “Diabetes Mellitus Type 2”). The graph made from the citation links has 438 connected components and the smallest positive eigenvalue of the augmented Normalized Laplacian is approximately µ̃ = 9.48 × 10−3. Therefore, the upper bound of Theorem 2 is λ−1 = (1− µ̃)−1 ≈ 1 + 9.57× 10−3. 7875 out of 19717 samples are labelled as “Diabetes Mellitus Type 2”, which is the largest proportion. Therefore, the chance rate is 7875/19717 = 39.9%.

heading: H.3.2 NOISY CITATION NETWORKS
text: We summarize the properties of noisy citation networks in Table 2.
12https://github.com/tkipf/gcn
We created two datasets from the Cora dataset: Noisy Cora 2500 and Noisy Cora 5000. Noisy Cora 2500 is made from the Cora dataset by uniformly randomly adding 2500 edges, respectively. Since some random edges are overlapped with existing edges, the number of newly-added edges is 2495 in total. We only changed the underlying graph from the Cora dataset and did not change word occurences (feature vectors) and genres (labels). The underlying graph of the Noisy Cora dataset has two connected components and the smallest positive eigenvalue is µ̃ ≈ 9.62 × 10−2. Therefore, the threshold of the maximum singular values of in Theorem 2 has been increased to λ−1 = (1− µ̃)−1 ≈ 1.11. Similarly, Noisy Cora 5000 was made by adding 5000 edges uniformaly randomly. The number of newly added edges is 4988 and the graph is connected (i.e., it has only 1 connected component). µ̃ and λ are µ̃ ≈ 1.32× 10−1 and λ = (1− µ̃)−1 ≈ 1.15, respectively. We made the noisy version of CiteSeer (Noisy CiteSeer) and PubMed (Noisy PubMed), in the similar way, by adding 5000 and 25000 edges uniformly randomly to the datasets. Since some random edges were overlapped with existing edges, 4991 and 24993 edges are newly added, respectively. This manipulation reduced the number of connected component of the graph to 3. µ̃ is approximately 1.11 × 10−1 (Noisy CiteSeer) and 1.43 × 10−1 (Noisy PubMed) and λ−1 = (1 − µ̃)−1 is approximately 1.13 (Noisy CiteSeer) and 1.17 (Noisy PubMed), respectively. Figure 5 (right) shows the spectral distribution of the augmented normalized Laplacian For comparison, we show in Figure 5 (left) the spectral distribution of the normalized Laplacian for these datasets13.

heading: H.3.3 MODEL ARCHITECTURE
text: We used a GCN consisting of a single node embedding layer, one to nine graph convolution layers, and a readout operation (Gilmer et al., 2017), which is a linear transformation common to all nodes in our case. We applied softmax function to the output of GCN. The output dimension of GCN is same as the number of classes (i.e., seven for Noisy Cora 2500/5000, six for Noisy CiteSeer, and three for Noisy PubMed). We treated the number of units in each graph convolution layer as a hyperparameter. Optionally, we specified the maximum singular values s of graph convolution
13Due to computational resource problems, we cannot compute the spectral distributions for PubMed and Noisy Pubmed.
layers. The choice of s is either 0.5 (smaller than 1), s1 (in the interval {1 ≤ s < λ−1}), 3 and 10 (larger than λ−1). We used s1 = 1.05 for Noisy Cora 2500, Noisy CiteSeer, and Noisy PubMed, and s1 = 1.1 for Noisy Cora 5000 and Noisy CiteSeer so that s1 is not close to the edges of the the interval {1 ≤ s < λ−1}.

heading: H.3.4 PERFORMANCE EVALUATION PROCEDURE
text: We split all nodes in a graph (either Noisy Cora 2500/5000 or Noisy CiteSeer) into training, validation, and test sets. Data split is the same as the one done by Kipf & Welling (2017). This is a transductive learning (Pan & Yang, 2010) setting because we can use node properties of the validation and test data during training. We trained the model three times for each choice of hyperparemeters using the training set and defined the objective function as the average accuracy on the validation set. We chose the combination of hyperparameters that achieves the best value of objective function. We evaluate the accuracy of the test dataset three times using the chosen combination of hyperparameters and computed their average and the standard deviation.

heading: H.3.5 TRAINING
text: At initialization, we sampled parameters from the i.i.d. Gaussian distribution. If the scale of maximum singular values s was specified, we subsequently scaled weight matrices of graph convolution layers so that their maximum singular values were normalized to s. The loss function was defined as the sum of the cross entropy loss for all training nodes. We train the model using the one of gradient-based optimization methods described in Table 3.

heading: H.3.6 HYPERPRAMETERS
text: Table 3 shows the set of hyperparameters from which we chose. Since we compute the representations of all nodes at once at each iteration, each epoch consists of 1 iteration. We employ Tree-structured Parzen Estimator (Bergstra et al., 2011) for hyperparameter optimization.
H.3.7 IMPLEMENTATION
We used Chainer Chemistry14, which is an extension library for the deep learning framework Chainer (Tokui et al., 2015; 2019), to implement GCNs and Optuna (Akiba et al., 2019) for hyperparameter tuning. We conducted experiments in a signel machine which has 2 Intel(R) Xeon(R) Gold 6136 CPU@3.00GHz (24 cores), 192 GB memory (DDR4), and 3 GPGPUs (NVIDIA Tesla V100). Our implementation achieved 68.1% with Dropout (Srivastava et al., 2014) (2 graph convolution layers) and 64.2% without Dropout (1 graph convolution layer) on the test dataset. These are slightly worse than the accuracy reported in Kipf & Welling (2017), but are still comparable with it.

heading: H.4 EXPERIMENT OF SECTION 6.4
text: The experiment settings are almost same as the experiment in Section 6.3. The only difference is that we did not train the node embedding layer, which we put before convolution layers of a GCN, while we did in Section 6.3. This is because we wanted to see the the effect of convolution operations
14https://github.com/pfnet-research/chainer-chemistry
on the perpendicular component of signals, while we interested in the prediction accuracy in real training settings in the previous experiment.

heading: I ADDITIONAL EXPERIMENT RESULTS

heading: I.1 EXPERIMENT OF SECTION 6.1
text: We show the vector field V for various W in Figure 6 (Case 1) and Figure 7 (Case 2). Parameters other than W are same as experiments in Section 6.1 (detail values are available in Appendix H.1).

heading: I.2 EXPERIMENT OF SECTION 6.2
text: Figure 8 shows the relative log distance of signals and their upper bound for various edge probability p and the maximum singular value s. Note that we generate a new graph for each configuration of (p, s). Therefore, different configurations may have different graphs and hence different λ even they have a same edge probability p in common.

heading: I.3 EXPERIMENT OF SECTION 6.3

heading: I.3.1 PREDICTIVE ACCURACY
text: Figure 9 shows the comparison of predictive performance in terms the maximum singular value and layer size when the dataset is Noisy Cora 5000 (left) and Noisy Citeseer (right), respectively. Concrete values are available in Table 4.

heading: I.3.2 TRANSITION OF MAXIMUM SINGULAR VALUES
text: Figure 10 – 13 show the transition of weight of graph convolution layers during training when the dataset is Noisy Cora 2500, Noisy Cora 5000, and Noisy CiteSeer, respectively. We note that the result of 3-layered GCN from the Noisy Cora 2500 is identical to Figure 3 (right) of the main article.

heading: I.4 EXPERIMENT OF SECTION 6.4
text: Figure 14 shows the logarithm of relative perpendicular component and prediction accuracy on Noisy Cora, Noisy CiteSeer, and Noisy PubMed datasets. We use Pearson R as a correlation coefficient. If GCNs have only one layer, it has more large relative perpendicular components (corresponding to right points in the figures) than GCNs which have other number of layers. The correlation between the logarithm of relative perpendicular components and prediction accuracies are 0.827(p = 6.890 × 10−6) for Noisy Cora, 0.524(p = 1.771 × 10−2) for Noisy CiteSeer, and 0.679(p = 1.002 × 10−3) for Noisy PubMed, if we treat the one-layer case as outliers and remove them.
dM(X(0)) and
orange dotted lines are upper bound y(l) := l log(sλ), where X(0) is the input signal and X(l) is the output of the l-th layer. Best view in color.

Review summaries:
In the paper, the authors carry out theoretically analysis on the expressive power for GNN. The analysis focused on the limiting case when the depth of layers goes to infinite. The authors prove that if the weights of the GNN satisfy certain condition based on the graph Laplacian, then the transformed features contain only degree and connected component information. Then the authors study the G_{np} random graph as a special case. Finally, empirical experiments are carried out to corroborates the theoretical results.

The paper studies why graph NNs lose the expressive power as additional layers are added. A dynamical system perspective is adopted and used to show that under certain conditions on the weights, the expressiveness of the network deteriorates. This is since the network's output eventually only carries information about superficial graph properties for distinguishing nodes, and nothing else. A case study of Erdos-Renyi graph is provided, showing that for dense graphs information loss indeed occurs. Guidelines to try and deal with this in practice were devised and empirically examined.