Title: GRAPH NEURAL NETWORKS EXPONENTIALLY LOSE EXPRESSIVE POWER FOR NODE CLASSIFICATION

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 data 1 .
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 Veličković 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

1. INTRODUCTION:

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.

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 Erdős -Rényi graph G N,p , which has N nodes and each edge appears independently with probability p, for an example. We prove that if log N pN = o(1) as a function of N , any GCN whose weights have maximum singular values at most C N p 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 Erdős -Ré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.

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 (JK-Net) 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.

3. PROBLEM SETTING AND MAIN RESULT:

3.1 NOTATION Let N + be the set of positive integers. For N ∈ N + , we denote [N ] := {1, . . . , N }. For a vector v ∈ R N , we write v ≥ 0 if and only if

v n ≥ 0 for all n ∈ [N ]. Similarly, for a matrix X ∈ R N ×C , we write X ≥ 0 if and only if X nc ≥ 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 ∈ R N and X, Y := tr(X T Y ) for X, Y ∈ R N ×C . 1 P equals to 1 if the proposition P is true else 0. For vectors v ∈ R N and w ∈ R C , v ⊗ w ∈ R N ×C denotes the Kronecker product of v and w defined by (v ⊗ w) nc := v n w c . For X ∈ R N ×C , X F := X, X 1 2 denotes the Frobenius norm of X. For a vector v ∈ R N , diag(v) := (v n δ nm ) n,m∈[N ] ∈ R N ×N

denotes the diagonalization of v. I N ∈ R N ×N denotes the identity matrix of size N . For a linear operator P : R N → R M and a subset V ⊂ R N , we denote the restriction of P to V by P | V .

3.2. DYNAMICAL SYSTEM:

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, H l ∈ N + (l ∈ N + ), let P ∈ R N ×N be a symmetric matrix and W lh ∈ R C×C for l ∈ N + and h ∈ [H l ]. We define f l : R N ×C → R N ×C by f l (X) := MLP l (P X). Here, MLP l : R N ×C → R N ×C is the l-th multi-layer perceptron common to all nodes (Xu et al., 2019) and is defined by

MLP l (X) := σ(• • • σ(σ(X)W l1 )W l2 • • • W lH l ), where σ : R N ×C → R N ×C

is an element-wise ReLU function (Krizhevsky et al., 2012) defined by σ(X) nc := max(X nc , 0)

for n ∈ [N ], c ∈ [C].

We consider the dynamics X (l+1) := f l (X (l) ) with some initial value X (0) ∈ R N ×C . We are interested in the asymptotic behavior of X (l) as l → ∞.

For M ≤ N , let U be a M -dimensional subspace of R N . 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 (e m ) m∈[M ] that consists of non-negative vectors. Assumption 2. U is invariant under P , i.e., if u ∈ U , then P u ∈ U .

We endow R N with the ordinal inner product and denote the orthogonal complement of U by

U ⊥ := {u ∈ R N | 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.

3.3. MAIN RESULT:

We define the subspace

M of R N ×C by M := U ⊗ R C = { M m=1 e m ⊗ w m | w m ∈ R C } where (e m ) m∈[M ] is the orthonormal basis of U appeared in Assumption 1. For X ∈ R N ×C , we denote the distance between X and M by d M (X) := inf{ X -Y F | Y ∈ M}.

We denote the maximum singular value of W lh by s lh and set s l := H l h=1 s lh . With these preparations, we introduce the main theorem of the paper. Theorem 1. Under Assumptions 1 and 2, we have

d M (f l (X)) ≤ s l λd M (X) for any X ∈ R N ×C .

The proof key is that the non-linear operation σ decreases the distance d M , that is, d M (σ(X)) ≤ d M (X). We use the non-negativity of e m 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 d M (X) = 0, this theorem implies that M is invariant under f l . 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 f l for any l ∈ N + , that is, if X ∈ M, then we have f l (X) ∈ M. Corollary 2. Let s := sup l∈N+ s l . We have d M (X (l) ) = O((sλ) l ). In particular, if sλ < 1, then X l exponentially approaches M 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 λ := max n∈[N -M ] |λ n |.

We denote U by the eigenspace associated with λ N and assume that U satisfies Assumption 1. Then, we have d M (X (l+1) ) ≤ s l λd M (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 W l = 1 for all l ∈ N + . This is essentially the direct consequence of Perron -Frobenius' theorem (see e.g., Meyer (2000) ). See Appendix F.

4. APPLICATION TO GCN:

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 in G by N = |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 ] ∈ R N ×N be the adjacency matrix and D := diag(deg(i) i∈[N ] ) ∈ R N ×N be the degree matrix of G where deg(i) := |{j ∈ V | (i, j) ∈ E}| is the degree of node i. Let Ã := A + I N , D := D + I N 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 ∆ := I N -D-1 2 Ã D-12 and set P := I N -∆. Let L, C ∈ N + be the layer and channel sizes, respectively. For weights

W l ∈ R C×C (l ∈ [L]), we define a GCN 2 associated with G by f = f L • • • • • f 1 where f l : R N ×C → R N ×C

is defined by f l (X) := σ(P XW l ). 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 = V 1 • • • V M be

the decomposition of the node set V into connected components. We denote an indicator vector of the m-th connected component by

u m := (1 {n∈Vm} ) n∈[N ] ∈ R N .

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 λ := max n=1,...,N -M |λ n | < 1. Further, e m := D 1 2 u m 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 d M (X (l) ) ≤ (sλ) l d M (X (0) ). In particular, d M (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.

5. ASYMPTOTIC BEHAVIOR OF GCN ON ERD ŐS -R ÉNYI GRAPH:

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 Erdős -Rényi graph G N,p (Erdös & Ré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 λ := min m=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. For an Erdős -Rényi graph G N,p , Chung & Radcliffe (2011) showed that when log N N p = 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 Erdős-Ré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 Erdős -Rényi graph. See Appendix D for the complete proof.

Theorem 3. Consider a GCN on the Erdős-Rényi graph G N,p such that log N N p = 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 < s 0 := 1 7 N p-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 = log N/N . Then, each node has the order of N p ≈ 4.3 adjacent nodes.

Under the condition of Theorem 3, the upper bound s 0 → ∞ 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, s 0 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).

6. EXPERIMENT 6.1 SYNTHESIS DATA: ONE-STEP TRANSITION:

We numerically investigate how the transition f (X) := σ(P XW ) 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 ∈ R 2 

(l) = log(d M (X (l) )/d M (X (0)

)) and 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. 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 = [e 1 e 2 ] such that the signs of e 1 and e 2 differ (Case 2), which violates Assumption 1. We see that M is not invariant under f and f does not uniformly decrease the distance from M. Therefore, we cannot remove the non-negativity assumption from Theorem 1.

6.2. SYNTHESIS DATA: DISTANCE TO INVARIANT SPACE:

We evaluate the distance to the invariant space M using synthesis data. We randomly generate an Erdős -Rényi graph, a GCN on it, and an input signal X (0) . We compute the distance between the lth intermediate output X (l) and the invariant space M for various edge probability p and maximum singular value s. 

(l) = log d M (X (l) ) d M (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 d M in such a case for future research.

6.3. REAL DATA: EFFECT OF MAXIMUM SINGULAR VALUES ON PERFORMANCE:

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. We can decompose the output X of a model as X = X 0 + X 1 (X 0 ∈ M, X 1 ∈ M ⊥ ). According to the theory, X 0 has limited information for node classification. We hypothesize that the model emphasizes the perpendicular component X 1 to perform good predictions. To quantitatively evaluate it, we define the relative magnitude of the perpendicular component of the output X by t(X) := X 1 /X 0 . 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 above5 .

6.4. REAL DATA: EFFECT OF SIGNAL COMPONENT PERPENDICULAR TO INVARIANT SPACE:



7. DISCUSSION:

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, Coja-Oghlan (2007) theoretically proved that if the expected average degree of G N,p is bounded, the smallest positive eigenvalue of the normalized Laplacian of G N,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 Graph-SAGE (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 space M. Our proposed weight scaling mechanism takes this strategy. Recently, Zhao & Akoglu (2020) has pro-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.

8. CONCLUSION:

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 Erdős -Ré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) 

A PROOF OF THEOREM 1:

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 λ = sup n∈[N -M ] |λ n | and s lh is the maximum singular value of W lh Lemma 1. For any X ∈ R N ×C , we have d . Then, we obtain the desired inequality as follows:

M (P X) ≤ λd M (X). Lemma 2. For any X ∈ R N ×C , we have d M (XW lh ) ≤ s lh d M (X). Lemma 3. For any X ∈ R N ×C , we have d M (σ(X)) ≤ d M (X).

d 2 M (P X) = N m=M +1 λ m w m 2 ≤ λ 2 N m=M +1 w m 2 ≤ λ 2 N m=M +1 w m 2 = λ 2 d 2 M (X).

Proof of Lemma 2. Using the same decomposition of X as the proof in Lemma 1, we have

XW lh = N m=1 e m ⊗ (W lh w m ) = M m=1 e m ⊗ (W lh w m ) + N m=M +1 e m ⊗ (W lh w m ).

Therefore, we have

d 2 M (XW lh ) = N m=M +1 W lh w m 2 ≤ s 2 lh N m=M +1 w m 2 = s 2 lh d 2 M (X).

Proof of Lemma 3. We choose (e m ) m=N -M +1,...,N as in the proof of Lemma 1. We denote X = (X nc ) n∈ 

(X) = N n=M +1 C c=1 a nc e c 2 ,

which is further transformed as

d 2 M (X) = N n=M +1 C c=1 a nc e c 2 = N n=M +1 C c=1 a 2 nc = C c=1 N n=1 a 2 nc - M n=1 a 2 nc = C c=1 X •c 2 - M n=1 X •c , e n 2 ,

where X •c is the c-th column vector of X. Similarly, we have

d 2 M (σ(X)) = C c=1 X + •c 2 - M n=1 X + •c , e n 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 ∈ R N and v 1 , . . . , v M ∈ R N be orthonormal vectors (i.e., v m , v n = δ mn ) satis- fying v m ≥ 0 for all m ∈ [M ]. Then, we have x 2 - M m=1 x, v m 2 ≥ x + 2 - M m=1 x + , v m 2 where x + := max(x, 0) for x ∈ R.

Proof. The value y 2 -M m=1 y, u m 2 is invariant under simultaneous coordinate permutation of y and u m 's. Therefore, we can assume without loss of generality that the coordinate of x are sorted:

EQUATION

When L = 0, the sum in the right hand side is treated as 0. On the other hand, writing as

EQUATION

Let

I m := {n ∈ [N ] | v nm > 0} be the support of v m for m ∈ [M ]. We note that if m = m ∈ [M ], we have I m ∩ I m = ∅ since if there existed n ∈ I m ∩ I m , we have 0 = v m , v m ≥ v nm v nm > 0,

which is contradictory. Therefore,

N m=1 L n=1 x n v nm 2 = N m=1   n∈Im∩[L] x n v nm   2 ≤ N m=1   n∈Im∩[L] x 2 n     n∈Im∩[L] v 2 nm   (∵ Cauchy-Schwarz inequality) ≤ N m=1   n∈Im∩[L] x 2 n   (∵ v m 2 = 1) ≤ L n=1

x 2 n .

(3)

We used the fact that I m 's are disjoint and v nm = 0 if n ∈ ∪ m I m in the first equality above. Further, we have x n x l v nm v lm ≤ 0 for 1 ≤ n ≤ L and L + 1 ≤ l ≤ N by the definition of L and non-negativity of v m . By combining (1), (2), and (3), we have

M m=1 x, v m 2 -x + , v m 2 ≤ L n=1 x 2 n = x 2 -x + 2 .

Proof of Theorem 1. By Lemma 1, 2, and 3, we have

d M (f l (X)) = d M (σ(• • • σ(σ( H times P X)W l1 )W l2 • • • W lH l )) ≤ d M (σ(• • • σ(σ( H-1 times P X)W l1 )W l2 • • • )W lH l ) ≤ s lH l -1 d M (σ(• • • σ(σ( H-1 times P X)W l1 )W l2 • • • )W lH l -1 )) • • • ≤ H l h=1 s lh d M (P X) ≤ s l d M (P X) ≤ s l λd M (X).

B PROOF OF PROPOSITION 1:

Proof. Let μ1 ≤ • • • ≤ μN be the eigenvalue of the augmented normalized Laplacian ∆, sorted in ascending order. Since P = I N -∆, 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 components6 . 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 a ij x i √ d i + 1 - x j d j + 1 2 for any x = [x 1 • • • x N ] ∈ R N .

Therefore, ∆ is positive semi-definite and hence μ1 ≥ 0.

Suppose temporally that G is connected. If x ∈ R N is an eigenvector associated to 0, then, by the aformentioned calculation, xi must be same for all pairs (i, j) such that a ij > 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 2 1. Now, suppose G has M connected components V 1 , . . . , V M . Let ∆m be the augmented normalized Laplacians corresponding to each connected component V m for m ∈ [M ]. By the aformentioned discussion, ∆m has the eigenvalue 0 with multiplicity 1. Since ∆ is the direct sum of ∆ m s, the eigenvalue of ∆ is the union of those for ∆m 's. Therefore, ∆ has the eigenvalue 0 with multiplicity M and e m = D 1 2 1 m 's are the orthogonal basis of the eigenspace.

Finally, we prove μN < 2. Let µ N be the largest eigenvalue of the normalized Laplacian

∆ = D -1 2 (D -A)D -1 2 , where D -1 2 ∈ R N ×N is the diagonal matrix defined by D -1 2 ii = deg(i) -1 2 (if deg(i) = 0) 0 (if deg(i) = 0) . Note that D -1 2 D 1 2 nor D 1 2 D -1 2

are not equal to the identity matrix I N in general. However, we have

EQUATION

where L = D -A is the (unnormalized) Laplacian. Therefore, we have

μN = max x =0 x ∆x x = max x =0 x D-1 2 L D-1 2 x x = max x =0 x D-1 2 D 1 2 D -1 2 LD -1 2 D 1 2 D-1 2 x x (∵ (4)) = max x =0 (D 1 2 D-1 2 x) ∆(D 1 2 D-1 2 x) x = max x =0 D 1 2 D-1 2 x =0 (D 1 2 D-1 2 x) ∆(D 1 2 D-1 2 x) x = max x =0 D 1 2 D-1 2 x =0 (D 1 2 D-1 2 x) ∆(D 1 2 D-1 2 x) D 1 2 D-1 2 x D 1 2 D-1 2 x x ≤ max x =0 D 1 2 D-1 2 x =0 (D 1 2 D-1 2 x) ∆(D 1 2 D-1 2 x) D 1 2 D-1 2 x max x =0 D 1 2 D-1 2 x =0 D 1 2 D-1 2 x x ≤ max y =0 y ∆y y max x =0 D 1 2 D-1 2 x x = µ N max n∈[N ] deg(i) deg(i) + 1 1 2 ≤ µ N . Therefore, we have μN ≤ µ N 7 . Since max i∈[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.

C COUNTEREXAMPLE OF PREVIOUS STUDY ON OVER-SMOOTHING FOR NON-LINEAR GNNS:

We restate Theorem 1 of the preprint (version2) of Luan et al. (2019) 8 . Let G 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/2 Ã D-1/2 ∈ R N ×N be the augmented normalized Laplacian of G. Let F ∈ N + and

W n ∈ R F ×F be the weight of the n-th layer for n ∈ N + . For the input X ∈ R N ×F , we define the output Y n ∈ R N ×F of the n-th layer of a GCN by

Y n = σ(L • • • σ(LXW 0 ) • • • W n )

where σ is the ReLU function. We assume the input X is drawn from a continuous distribution on R N ×F . Then, the theorem claims that we have lim n→∞ rank(Y n ) = 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 1 1 1 1 0 1 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 W n = I 3 (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 R 4×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 Y n = L n X 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(Y n ) = 3 ⇐⇒ rank(P Λ n P X) = 3 ⇐⇒ X ∈ {P -1 [B v] | B ∈ R 3×3 is invertible, v ∈ R 3 }.

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(Y n ) = 3 for all n ∈ N}) > 0. Therefore, we have lim n→∞ rank(Y n ) = 3 with a non-zero probability.

D PROOF OF THEOREM 3:

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 Erdős-Ré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 Erdős-Rényi graphs. Let N ∈ N + , and P = (p ij ) i,j∈[N ] be a non-negative symmetric matrix (meaning that p ij ≥ 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 p ij . Let A and D be the adjacency and the degree matrices of G, respectively (that is, A ij ∼ Ber(p ij ), i.i.d.). Define the expected node degree of node i by t i := 

, if δ + 1 > k(ε) log N , we have max n∈[N ] λ n ( ∆) -λ n ( ∆) ≤ 4 3 log(4N/ε) δ + 1

with probability at least 1 -ε.

Proof. By Weyl's theorem, we have max n∈

[N ] λ n ( ∆) -λ n ( ∆) ≤ ∆ -∆ . Therefore, it is enough to bound ∆ -∆ . Let C := I N -D-1 2 Ã D.

By the triangular inequality, we have ∆ -∆ ≤ ∆ -C + C -∆ . We will bound these terms respectively.

First, we bound C-∆ . Direct calculation shows C-∆ = -D-1 2 (A-P ) D-1 2 . Let E ij ∈ R N ×N be a matrix defined by

(E ij ) kl = 1 if (i = k and i = l) or (i = l and j = k), 0 otherwise.

We define the random variable Y ij by

Y ij := A ij -p ij √ t i + 1 t j + 1 E ij .

Then, Y ij 's are independent and we have

C -∆ = N i,j=1 Y ij . To apply Theorem 5 of Chung & Radcliffe (2011) to Y ij 's, we bound Y ij -E[Y ij ] and N i,j=1 E[Y 2 ij ] . First, we have Y ij -E[Y ij ] = Y ij ≤ E ij √ t i + 1 t j + 1 ≤ (δ + 1) -1 . Since E[Y 2 ij ] = p ij -p 2 ij (t i + 1)(t j + 1) E ii + E jj (if i = j), E ii (if i = j),

we have

N i,j=1 E[Y 2 ij ] = N i,j=1 p ij -p 2 ij (t i + 1)(t j + 1) E ii = max i∈[N ]   N j=1 p ij -p 2 ij (t i + 1)(t j + 1)   ≤ max i∈[N ]   N j=1 p ij (t i + 1)(t j + 1)   ≤ (δ + 1) -1 .

By letting a ← 3 log(4N/ε) δ+1 , M ← (δ + 1) -1 , v 2 ← (δ + 1) -1 and applying Theorem 5 of Chung & Radcliffe (2011) , we have

Pr( C -∆ > a) ≤ 2N exp - a 2 2(δ + 1) -1 + 2(δ + 1) -1 a/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

EQUATION

)

Next, we bound ∆ -C . First, since a < 1, by Chernoff bound (see, e.g. Angluin & Valiant (1979) ; Hagerup & Rüb (1990) )), we have

Pr(|d i -t i | > a(t i + 1)) ≤ 2 exp - a 2 (t i + 1) 3 ≤ 2 exp - a 2 (δ + 1) 3 = ε 2N .

Therefore, if |d i -t i | ≤ a(t i + 1), then we have

d i + 1 t i + 1 -1 ≤ d i + 1 t i + 1 -1 (∵ | √ x -1| ≤ |x -1| for x ≥ 0) = d i -t i t i + 1 ≤ a.

Therefore, by union bound, we have

D-1 2 D 1 2 -I N = max i∈[N ] d i + 1 t i + 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

EQUATION

From ( 5) and ( 6), we have

∆ -∆ ≤ ∆ -C + C -∆ ≤ a + a(a + 1) + a ≤ a 2 + 3a ≤ 4a (∵ a < 1)

with probability at least 1 -ε by union bound.

Let N ∈ N + and p > 0. In the case of the Erdős-Rényi graph G N,p , we should set P = p(J N -I N ) where J N ∈ R N ×N are the all-one matrix. Then, we have

Ā = pJ N + (1 -p)I N , D = (N p - p + 1)I N , and ∆ = p N p-p+1 (N I N -J N ).

Since the eigenvalue of J N is N (with multiplicity 1) and 0 (with multiplicity N -1), the eigenvalue of ∆ is 0 (with multiplicity 1) and N p N p-p+1 (with multiplicity N -1). For G N,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 Erdős-Rényi graph G N,p . For any

ε > 0, if N p-p+1 log N > k(ε) := 3(1 + log(4/ε)), then, with probability at least 1 -ε, we have max i=2,...,N λ i ( ∆) - N p N p -p + 1 ≤ 4 3 log(4N/ε) N p -p + 1 .

Corollary 4. Consider GCN on G N,p . Let W l be the weight of the l-th layer of GCN and s l be the maximum singular value of W l for l ∈ N + . Set s := sup l∈N+ . Let ε > 0. We define k(ε

) := 3(1 + log(4/ε)) and l(N, p, ε) = 1-p N p-p+1 + 4 3 log(4N/ε) N p-p+1 . If N p-p+1 log N

> k(ε) and s ≤ l(N, ε) -1 , then, GCN on G N,p satisfies the assumption of Theorem 2 with probability at least 1 -ε.

Proof of Theorem 3. Since log N N p = o(1), for fixed ε, we have

N p -p + 1 log N > N p log N > k(ε)

for sufficiently large N . Further, N p → ∞ as N → ∞ when log N N p = o(1). Therefore, we have

(1 -p) 2 N p -p + 1 ≤ 1 N p ≤ (7 -4 √ 3) 2 log 4N ε

for sufficiently large N . Hence.

1 -p N p -p + 1 ≤ (7 -4 √ 3) log(4N/ε) N p -p + 1 .

Therefore, we have l(N, p, ε)

≤ 7 log(4N/ε) N p-p+1 . Therefore, if s ≤ 1 7 N p-p+1

log(4N/ε) , then we have s ≤ l(N, p, ε) -1 .

E MISCELLANEOUS PROPOSITIONS E.1 INVARIANCE OF ORTHOGONAL COMPLEMENT SPACE:

Proposition 2. Let P ∈ R N ×N be a symmetric matrix, treated as a linear operator P :

R N → R N . If a subspace U ⊂ R N is invariant under P (i.e., if u ∈ U , then P u ∈ U ), then, U ⊥ is invariant under P , too.

Proof. For any u ∈ U ⊥ and v ∈ U , by symmetry of P , we have

P u, v = (P u) v = u P v = u P v = u, P v .

Since U is an invariant space of P , we have P v ∈ U . Hence, we have u, P v = 0 because u ∈ U ⊥ . We obtain P u ∈ U ⊥ by the definition of U ⊥ .

E.2 CONVERGENCE TO TRIVIAL FIXED POINT:

Let P ∈ R N ×N be a symmetric matrix, W l ∈ R C×C , s l be the maximum singular value of W l for l ∈ N + . We define f l : R N ×C → R N ×C by f l (X) := σ(P XW l ) where σ is the element-wise ReLU function. Proposition 3. Suppose further that the operator norm of P is no larger than λ, then we have f l (X) F ≤ s l λ X F for any l ∈ N + . In particular, let s := sup l∈N+ s l . If sλ < 1, then, X l 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

P XW l F ≤ λ XW l F ≤ s l λ X F . On the other hand, since σ(x) 2 ≤ x 2 for any x ∈ R, we have σ(X) F ≤ X F for any X ∈ R N ×C . Combining the two, we have f l (X) F ≤ P XW l F ≤ s l λ X F .

E.3 STRICTNESS OF THEOREM 1:

Theorem 1 implies that if sλ ≤ 1, then, one-step transition f l does not increase the distance to M. In this section, we first prove that this theorem is strict in the sense that, there exists a situation in which s l λ > 1 holds and the distance d M increases by one-step transition f l 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 : R N ×C → R N ×C by f (X) := σ(P XW ) where σ is the element-wise ReLU function.

Proposition 4. We have d

M (f (X)) > d M (X) for any X = [x 1 x 2 ] ∈ R 2 such that x 2 > 0.

Proof. By definition, we have d M (X) = |x 2 |. On the other hand, direct calculation shows that f l (X) = (W µX 1 ) + (W λX 2 ) + and d M (f l (X)) = (W λX 2 ) + where x + := max(x, 0) for

x ∈ R. Since W > λ -1 and x 2 > 0, we have d M (f l (X)) > d M (X).

Next, we prove the non-strictness of Theorem 1 in the sense that there exists a situation in which s l λ > 1 holds and the distance d M uniformly decreases by f l . 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 d M (f l (X)) < d M (X) for all X ∈ R 2 . Proof. First, note that e := 1 √ 2 [1 -1]

is the eigenvector of P associated to λ: P e = λe . For X = ae + be (a, b ¿ 0), the distance between X and M is d M (X) = |b|. On the other hand, by direct computation, we have

f (X) = σ(P XW ) =      0 λb √ 2 (if b ≥ 0), -λb √ 2 0 (if b < 0).

Therefore, the distance between f (X) and

M is d M (f (X)) = λ|b|/2. Since λ < 2, we have d M (f (X)) < d M (X) for any X ∈ R 2 .

We have shown that the non-negativity of e (Assumption 1) is not a redundant condition in Section 6.1.

F RELATION TO MARKOV PROCESS:

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 = (p ij ) i,j∈[N ] ∈ R N ×N such that P ≥ 0 and P 1 = 1 where 1 is the all-one vector. We interpret p ij as the transition probability from a state i to j. We associate P with a graph G P = (V P , E P ) by V P = [N ] and (i, j) ∈ E P if and only if p ij > 0. Since P is symmetric, we can regard G P as an undirected graph. We assume P is irreducible and aperiodic10 . 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 ∈ R N such that e ≥ 0, e = 1, and e is the eigenvector for the eivenvalue 1. Corollary 5. Let λ := max n=1,...,N -1 |λ n |(< 1) and M := {e ⊗ w | w ∈ R C }. If sλ < 1, then, for any initial valueX 1 , X l exponentially approaches M as l → ∞.

If we set C = 1 and W l = 1 for all l ∈ N + , then, we can inductively show that X l ≥ 0 for any l ≥ 2. Therefore, we can interpret X l as a measure on S. Suppose further that we take the initial value X 1 as X 1 ≥ 0 and X 1 1 = 1 so that we can interpret X 1 as a probability distribution on S. Then, we can inductively show that X l ≥ 0, X l 1 = 1 (i.e., X l is a probability distribution on S), and X l+1 = σ(P X l W l ) = P X l 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.

G GCN DEFINED BY NORMALIZED LAPLACIAN:

In Section 4, we defined P using the augmented normalized Laplacian ∆ by P = I N -∆. 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 ∆ := I N -∆. 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) ).

H DETAILS OF EXPERIMENT SETTINGS:

H.1 EXPERIMENT OF SECTION 6.1

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. We randomly generated an Erdős -Rényi graph G N,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) ∈ R N ×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). 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%.

H.3.2 NOISY CITATION NETWORKS:

We summarize the properties of noisy citation networks in Table 2 . 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 .

H.3.3 MODEL ARCHITECTURE:

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 layers. The choice of s is either 0.5 (smaller than 1), s 1 (in the interval {1 ≤ s < λ -1 }), 3 and 10 (larger than λ -1 ). We used s 1 = 1.05 for Noisy Cora 2500, Noisy CiteSeer, and Noisy and s 1 = 1.1 for Noisy Cora 5000 and Noisy CiteSeer so that s 1 is not close to the edges of the the interval {1 ≤ s < λ -1 }.

H.3.4 PERFORMANCE EVALUATION PROCEDURE:

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.

H.3.5 TRAINING:

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 .

H.3.6 HYPERPRAMETERS:

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.

H.4 EXPERIMENT OF SECTION 6.4

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 on the perpendicular component of signals, while we interested in the prediction accuracy in real training settings in the previous experiment.

I ADDITIONAL EXPERIMENT RESULTS:

I.1 EXPERIMENT OF SECTION 6.1

We show the vector field V for various W in Figure 6 (Case 1) and Figure 7 I.4 EXPERIMENT OF SECTION 6.4

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. 

:

Code is available at https://github.com/delta2323/gnn-asymptotics.

Following the original paper(Kipf & Welling, 2017), we use one-layer MLPs (i.e., H l = 1 for all l ∈ N+.). However, our result holds for the multi-layer case

Chung et al. (2004) and Coja-Oghlan (2007) proved similar theorems.

Since we consider the one-step transition only, we omit the subscript l from f l , X l , and W l .

We cannot conclude that large perpendicular components are essential for good performance, since the maximum singular value s is correlated to the accuracy, too.

The former statement is identical to Lemma 1 and latter one is the extension of Lemma 2 ofWu et al. (2019a).

Theorem 1 ofWu et al. (2019a) showed that this inequality strictly holds when G is simple and connected. We do not require this assumption.

https://arxiv.org/abs/1906.02174v2

Note that E[ ∆] = ∆ in general due to the dependence between Ã and D.

A 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.

This is the default initialization method for weight matrices in Chainer and Chainer Chemistry.

https://github.com/tkipf/gcn

Due to computational resource problems, we cannot compute the spectral distributions for PubMed and Noisy Pubmed.

https://github.com/pfnet-research/chainer-chemistry

Figures and Tables:

Figure 1: Visualization of vector field V (X) := f (X)-X induced by the one-step transition. Color maps indicate the absolute value |V (X)| at the point X. Dotted lines are the subspace M. Left: Case 1. Right: Case 2. Best view in color.

Figure2: The actual distances to the invariant space M and their upper bounds. Solid lines are the log relative distance defined by y(l) = log(d M (X (l) )/d M (X (0) )) and 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.

Figure 3: Node prediction results on Noisy Cora. Left: Effect of the maximum singular values on weights on model performance. The horizontal dotted line indicates the chance rate (30.2%). The error bar is the standard deviation of 3 trials. Right: Transition of maximum singular values during training. See Appendix I.3 for results using other datasets. Best view in color.

Figure 2 plots the logarithm of the relative distance y

Figure 4: log t(X) and prediction accuracy on Noisy Cora.

Proof of Lemma 1. Since P is a symmetric linear operator on U ⊥ , we can choose the orthonormal basis (e m ) m=M +1,...,N of U ⊥ consisting of the eigenvalue of P | U ⊥ . Let λ m be the eigenvalue of P to which e m 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 (e m ) m∈[N ] forms the orthonormal basis of R N , we can uniquely write X ∈ R N ×C as X = N m=1 e m ⊗ w m for some w m ∈ R C . Then, we have d 2 M (X) = N m=M +1 w m 2 where • is the 2-norm of a vector. On the other hand, we have (λ m w m ) Since U is invariant under P , for any m ∈ [M ], we can write P e m as a linear combination of e n (n ∈ [M ]). Therefore, we have d 2 M (P X) = N m=M +1 λ m w m 2

[N ],c∈[C] and e n = (e mn ) m∈[N ] , respectively. Let (e c ) c∈[C] be the standard basis of R C . Then, (e n ⊗ e c ) n∈[N ],c∈[C] is the orthonormal basis of R N ×C , endowed with the standard inner product as a Euclid space. Therefore, we can decompose X as X = N n=1 C c=1 a nc e n ⊗ e c where a nc = X, e n ⊗ e c = N m=1 X mc e mn . Then, we have d 2 M

j=1 p ij . Let Ã := A + I N , D := D + I N and define Ā := E[ Ã] = P + I N and D := E[ D] = diag(t 1 , . . . , t N ) + I N correspondingly. We define the augmented normalized Laplacian ∆ of G by ∆ := I N -D-1 2 Ã D-1 2 and its expected version by ∆ := I Na symmetric matrix X ∈ R N , we define its eigenvalues, sorted in ascending order by λ 1 (X) ≤ • • • ≤ λ N (X) and its operator norm by X = max n∈[N ] |λ n (X)|. Lemma 5 (Ref. Chung & Radcliffe (2011) Theorem 2). Let δ := min n∈[N ] t n be the minimum expected degree of G. Set k(ε) := 3(1 + log(4/ε)). Then, for any ε > 0

Figure 5: Spectral distribution of Laplacian for the citation network datasets. Left: normalized Laplacian. Right: augmented normalized Laplacian. Top: Cora and Noisy Cora (2500, 5000). Bottom: CiteSeer and Noisy CiteSeer.

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.

Figure 6: Vector field V for various W for Case 1. Top left: W = 0.5. Top right: W = 1.0. Middle left: W = 1.2 (same as Figure 1 in the main article). Middle right: W = 1.5. Bottom left: W = 2.0. Bottom right: W = 4.0. Best view in color.

1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1

Figure 7: Vector field V for various W for Case 2. Top left: W = 0.5. Top right: W = 1.0. Middle left: W = 1.2 (same as Figure 1 in the main article). Middle right: W = 1.5. Bottom left: W = 2.0. Bottom right: W = 4.0. Best view in color.

Figure 9: Effect of the maximum singular values of weights on predictive performance. Horizontal dotted lines indicate the chance rates (30.2% for Noisy Cora 5000, 21.2% for Noisy CiteSeer, and 39.9% for Noisy PubMed). The error bar is the standard deviation of 3 trials. Left: Noisy Cora 5000. Right: Noisy CiteSeer. Bottom: Noisy Pubmed. Best view in color.

Figure 10: Transition of maximum singular values of GCN during training using Noisy Cora 2500. Top left: 1 layer. Top right: 5 layers. Bottom left: 7 layers. Bottom right: 9 layers.

Figure 11: Transition of maximum singular values of GCN during training using Noisy Cora 5000. Top left: 1 layer. Top right: 3 layers. Middle left: 5 layers. Middle right: 7 layers. Bottom: 9 layers.

Figure 12: Transition of maximum singular values of GCN during training using Noisy CiteSeer. Top left: 1 layer. Top right: 3 layers. Middle left: 5 layers. Middle right: 7 layers. Bottom: 9 layers.

of graph NNs via spectral and random graph theories. Seiya Tokui, Ryosuke Okuta, Takuya Akiba, Yusuke Niitani, Toru Ogawa, Shunta Saito, Shuji Suzuki, Kota Uenishi, Brian Vogel, and Hiroyuki Yamazaki Vincent. Chainer: A deep learning framework for accelerating the research cycle. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pp. 2002-2011. ACM, 2019. Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In International Conference on Learning Representations, 2018. Saurabh Verma and Zhi-Li Zhang. Stability and generalization of graph convolutional neural networks. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, 2019. Boris Weisfeiler and Lehman A.A. A reduction of a graph to a canonical form and an algebra arising during this reduction. Nauchno-Technicheskaya Informatsia, 2(9):12-16, 1968. Felix Wu, Tianyi Zhang, Amauri Holanda de Souza Jr, Christopher Fifty, Tao Yu, and Kilian Q Weinberger. Simplifying graph convolutional networks. arXiv preprint arXiv:1902.07153, 2019a. Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and Philip S Yu. A comprehensive survey on graph neural networks. arXiv preprint arXiv:1901.00596, 2019b. Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanie Jegelka. Representation learning on graphs with jumping knowledge networks. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 5453-5462. PMLR, 2018. Pan Zhou and Jiashi Feng. Understanding generalization and optimization performance of deep CNNs. In Proceedings of the 35th International Conference on Machine Learning, volume 80 of Proceedings of Machine Learning Research, pp. 5960-5969. PMLR, 2018.

Table content: <table><tr><td>Ding-Xuan Zhou.</td><td>Universality of deep convolutional neural networks.</td><td>arXiv preprint</td></tr><tr><td colspan="2">arXiv:1805.10769, 2018.</td><td/></tr></table>

Dataset specifications. The threshold λ -1 in the table indicates the upper bound of Corollary 2. #Node #Edge #Class Chance Rate Threshold λ -1

Table content: <table><tr><td>Cora</td><td>2708</td><td>5429</td><td>6</td><td>30.2%</td><td>1 + 3.62 × 10 -3</td></tr><tr><td colspan="2">CiteSeer 3312</td><td>4732</td><td>7</td><td>21.1%</td><td>1 + 1.25 × 10 -3</td></tr><tr><td colspan="4">PubMed 19717 44338 3</td><td>39.9%</td><td>1 + 9.57 × 10 -3</td></tr></table>

Dataset specifications for noisy citation networks. The threshold λ -1 in the table indicates the upper bound of Corollary 2.Original Dataset #Edge Added Threshold λ -1

Table content: <table><tr><td colspan="2">Noisy Cora 2500 Cora</td><td>2495</td><td>1.11</td></tr><tr><td colspan="2">Noisy Cora 5000 Cora</td><td>4988</td><td>1.15</td></tr><tr><td>Noisy CiteSeer</td><td>CiteSeer</td><td>4991</td><td>1.13</td></tr><tr><td>Noisy PubMed</td><td>PubMed</td><td>24993</td><td>1.17</td></tr><tr><td colspan="2">H.3 EXPERIMENT OF SECTION 6.3</td><td/><td/></tr><tr><td>H.3.1 DATASET</td><td/><td/><td/></tr><tr><td>We used the Cora (McCallum</td><td/><td/><td/></tr></table>

Hyperparameters of the experiment in Section 6.3. X ∼ LogUnif[10 a , 10 b ] denotes the random variable log 10 X obeys the uniform distribution over[a, b]. "Learning rate" corresponds to α when "Optimization algorithm" is Adam(Kingma & Ba, 2015).

Table content: <table><tr><td>Name</td><td>Value</td></tr><tr><td>Unit size</td><td>{10, 20, . . . , 500}</td></tr><tr><td>Epoch</td><td>{10, 20, . . . , 100}</td></tr><tr><td colspan="2">Optimization algorithm {SGD, MomentumSGD, Adam}</td></tr><tr><td>Learning rate</td><td>LogUnif[10 -5 , 10 -2 ]</td></tr></table>

(Case 2). Parameters other than W are same as experiments in Section 6.1 (detail values are available in Appendix H.1).Figure8shows 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. 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 Table4.

Table content: <table><tr><td>I.2 EXPERIMENT OF SECTION 6.2</td></tr><tr><td>I.3 EXPERIMENT OF SECTION 6.3</td></tr><tr><td>I.3.1 PREDICTIVE ACCURACY</td></tr><tr><td>Figure 9 shows</td></tr></table>

Comparison of performance in terms of maximum singular value of weights and layer size. "U" in the right most column indicates the accuracy of GCN without weight normalization. ± 0.101 0.429 ± 0.090 0.552 ± 0.014 0.632 ± 0.007 0.587 ± 0.008 3 0.273 ± 0.051 0.309 ± 0.017 0.580 ± 0.058 0.661 ± 0.003 0.494 ± 0.041 5 0.319 ± 0.000 0.267 ± 0.059 0.462 ± 0.065 0.602 ± 0.004 0.326 ± 0.029 7 0.261 ± 0.076 0.262 ± 0.080 0.407 ± 0.021 0.501 ± 0.017 0.279 ± 0.129 9 0.261 ± 0.080 0.319 ± 0.000 0.284 ± 0.109 0.443 ± 0.014 0.319 ± 0.000 ± 0.080 0.333 ± 0.099 0.557 ± 0.004 0.561 ± 0.019 0.555 ± 0.016 3 0.245 ± 0.066 0.247 ± 0.076 0.370 ± 0.041 0.587 ± 0.009 0.286 ± 0.066 5 0.274 ± 0.048 0.237 ± 0.070 0.257 ± 0.076 0.535 ± 0.031 0.319 ± 0.000 7 0.263 ± 0.080 0.297 ± 0.031 0.260 ± 0.074 0.339 ± 0.060 0.319 ± 0.000 9 0.262 ± 0.081 0.258 ± 0.064 0.262 ± 0.080 0.261 ± 0.082 0.318 ± 0.002 ± 0.018 0.467 ± 0.012 0.490 ± 0.016 0.494 ± 0.006 0.495 ± 0.009 3 0.438 ± 0.027 0.436 ± 0.010 0.450 ± 0.019 0.462 ± 0.007 0.417 ± 0.061 5 0.285 ± 0.008 0.371 ± 0.016 0.373 ± 0.011 0.425 ± 0.007 0.380 ± 0.024 7 0.213 ± 0.006 0.282 ± 0.011 0.309 ± 0.012 0.385 ± 0.007 0.308 ± 0.012 9 0.182 ± 0.005 0.242 ± 0.030 0.303 ± 0.021 0.325 ± 0.003 0.229 ± 0.033 ± 0.039 0.636 ± 0.006 0.641 ± 0.010 0.632 ± 0.002 0.631 ± 0.010 3 0.442 ± 0.027 0.426 ± 0.026 0.658 ± 0.004 0.661 ± 0.005 0.631 ± 0.013 5 0.431 ± 0.033 0.431 ± 0.034 0.561 ± 0.083 0.641 ± 0.004 0.424 ± 0.093 7 0.428 ± 0.032 0.443 ± 0.051 0.449 ± 0.035 0.619 ± 0.011 0.440 ± 0.041 9 0.413 ± 0.009 0.438 ± 0.039 0.539 ± 0.052 0.569 ± 0.042 0.473 ± 0.031

Table content: <table><tr><td/><td/><td>Noisy Cora 2500</td><td/><td/></tr><tr><td/><td/><td colspan="2">Maximum Singular Value</td><td/></tr><tr><td>Depth 1</td><td>1.05</td><td>3</td><td>10</td><td>U</td></tr><tr><td colspan="3">1 0.389 Noisy Cora 5000</td><td/><td/></tr><tr><td/><td/><td colspan="2">Maximum Singular Value</td><td/></tr><tr><td>Depth 1</td><td>1.1</td><td>3</td><td>10</td><td>U</td></tr><tr><td colspan="3">1 0.301 Noisy CiteSeer</td><td/><td/></tr><tr><td/><td/><td colspan="2">Maximum Singular Value</td><td/></tr><tr><td>Depth 0.5</td><td>1.1</td><td>3</td><td>10</td><td>U</td></tr><tr><td colspan="3">1 0.461 Noisy Pubmed</td><td/><td/></tr><tr><td/><td/><td colspan="2">Maximum Singular Value</td><td/></tr><tr><td>Depth 0.5</td><td>1.1</td><td>3</td><td>10</td><td>U</td></tr><tr><td>1 0.488</td><td/><td/><td/><td/></tr></table>

Peer review (iclr2020):
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. I enjoyed reading the paper, and found the results interesting and informative. The presentation is good and the derivations are clear. I do not have significant criticism, but would like address a few points (see below). In sum, the paper should be accepted in my opinion. Comments: 1) Understanding the behavior in sparse graphs seems like an important research avenue, without which the picture is incomplete. Perhaps adding some empirical results in that direction? 2) The method for transforming sparse graphs into dense ones in Section 6.3 is not convincing. The random addition of edges compromises the truth correlation profile of citations. What should we learn about the actual task from the experiments on the noisy version? 3) In relation to generalization error, mentioned in the discussion section. Could the double descent curves of https://arxiv.org/abs/1812.11118 reconcile the discussion therein? How would over-parameterization from the interpolation perspective and the results of the current work relate?