Revisiting "Over-smoothing" in Deep GCNs

1 Introduction

Graph neural networks (GNNs) are widely used in modeling real-world connections, like protein networks

ying2018hierarchical , social networks velivckovic2018deep , and co-author networks kipf2016semi . One could also construct similarity graphs by linking data points that are close in the feature space even when there is no explicit graph structure. There have been several successful GNN architectures: ChebyshevNet defferrard2016convolutional , GCN kipf2016semi , SGC wu2019simplifying , GAT velivckovic2017graph , GraphSAGE hamilton2017inductive and other subsequent variants tailored for practical applications gilmer2017neural ; liu2019hyperbolic ; klicpera2018predict .

Recently, researchers started to explore the fundamentals of GNNs, such as expressive power xu2018powerful ; oono2019graph ; loukas2019graph ; dehmamy2019understanding , and analyze their capacity and limitations. One of the frequently mentioned limitations is oversmoothing li2019deepgcns . In deep GCNs, over-smoothing means that after multi-layer graph convolution, the effect of Laplacian smoothing makes node representations become more and more similar, eventually becoming indistinguishable. This issue was first mentioned in li2018deeper and has been widely discussed since then, such as in JKNet xu2018representation , DenseGCN li2019deepgcns , DropEdge rong2019dropedge , and PairNorm zhao2019pairnorm . However, these discussions were mostly on the powering effect of convolution operator $A^{L}$ (where $A$ is the convolution operator, and $L$

is the number of layers). This essentially implies a GCN variant without an activation function (i.e., a simplified graph convolution, or SGC

wu2019simplifying ).

In this work, we instead argue that deep GCNs can learn anti-oversmoothing, while overfitting is the real cause of performance drop. This paper focuses on graph node classification and starts from the perspective of graph-based optimization (minimizing $L_{0} + γ L_{r e g}$ ) belkin2003laplacian ; yang2016revisiting , where $L_{0}$ is the (supervised) empirical loss and $L_{r e g}$ is an (unsupervised) graph regularizer, which encodes smoothness over the connected node pairs. We interpret a GCN as a two-step optimization problem: (i) forward propagation to minimize $L_{r e g}$ by viewing ${W^{(l)}}$ as constants and ${X^{(l)}}$ as parameters, and (ii) back propagation to minimize $L_{0}$ by updating ${W^{(l)}}$ . This work therefore derives GCN as:

STEP1: min {X^{(l)}} L_{r e g} ({X^{(l)}} ∣ {W^{(l)}}) and STEP2: min {W^{(l)}} L_{0} ({W^{(l)}})

Similarly, an SGC is interpreted as:

STEP1: min X L_{r e g} (X) and STEP2: min W L_{0} (W)

From this formulation, we show that a deep SGC indeed suffers from oversmoothing but deep GCNs will learn to prevent oversmoothing because (i) $L_{r e g}$ is conditioned on ${W^{(l)}}$ in GCN, and (ii) $L_{0}$ and $L_{r e g}$ are somehow contradicted. To prevent gradient vanishing/exploding, in this paper, we add skip connections he2016deep to all deep architectures by default. An illustration of deep GCNs, deep SGC and directly learning $L_{0} + γ L_{r e g}$ using DNN is shown in Fig. 1.

[width=5.5in]fig/GCNvsSGC.pdf

Figure 1: Comparison of Deep GCN, Deep SGC and DNN on Cora

As is mentioned above, the training of deep GCNs is a learning process of anti-oversmoothing, which is extremely slow in practice and sometimes may not converge. Based on the formulation, we further propose a mean-subtraction trick to accelerate the training of deep GCNs. Extensive experiments verify our theories and provide more insights about deep GCNs.

2 Background of Graph Transductive Learning

Graph representation learning aims at embedding the nodes into low-dimensional vectors, while simultaneously preserving both

graph topology structure and node feature information. Given a graph

G = (V, E)

, let

V = {v_{1}, v_{2}, \dots, v_{n}}

be the set of nodes, and let

Y

be a set of

m

possible classes. Assume that each node

v_{j}

is associated with a class label

y_{j} \in Y

. A graph could be represented by an adjacency matrix

A

with

A_{i j} = 1

when two nodes are connected

(v_{i}, v_{j}) \in E

. The degree matrix

D = d i a g (d_{1}, d_{2}, \dots, d_{n})

is diagonal where

d_{i} = \sum_{j} A_{i j}

. Let

X = {x_{1}, x_{2}, \dots, x_{n}}

denote the feature vectors for each node. Given a labelled set

T \subset V

, the goal of transductive learning on a graph is to transductively predict labels for the remaining unknown nodes

V ∖ T

. A well-studied solution category is to include graph regularizers tenenbaum2000global ; belkin2003laplacian ; zhu2002learning ; weston2012deep into the classification algorithm. Graph-convolution-based models kipf2016semi ; defferrard2016convolutional ; velivckovic2017graph ; hamilton2017inductive are a powerful learning approach in this space.

2.1 Graph-based Regularization

There is a rather general class of embedding algorithms that include graph regularizers. They could be described as: finding a mapping $f (\cdot)$ by minimizing the following two-fold loss:

L = L_{0} (f (X)) + γ L_{r e g} (f (X)),

(1)

where $f (X) = [f (x_{i})]_{i = 1}^{n}$ is the low-dimensional representation of nodes. The first term is the empirical risk on the labelled set $T$ . The second term is a graph regularizer over connected pairs, so as to make sure that a trivial solution is not reached.

The measurements on graphs are usually invariant to node permutations. A canonical way is to use Dirichlet energy belkin2002laplacian for the graph-base regularization,

L_{r e g} = \frac{1}{2} \sum i, j A_{i j} {∥ ∥ ∥ ∥ \frac{f (x_{i})}{\sqrt{d_{i}}} - \frac{f (x_{j})}{\sqrt{d_{j}}} ∥ ∥ ∥ ∥}^{2} = \frac{1}{2} Tr (f (X)^{⊤} Δ f (X)),

(2)

where $Δ = I - D^{- \frac{1}{2}} A D^{- \frac{1}{2}}$ is the normalized Laplacian operator, which induces a semi-norm on $f (\cdot)$ , penalizing the changes between adjacent vertices. Same normalized formulation could be found in chen2019deep ; ando2007learning ; smola2003kernels ; shaham2018spectralnet ; belkin2003laplacian , and some related literature also use the unnormalized version zhu2002learning ; von2007tutorial .

2.2 Graph Convolutional Network

GCNs are derived from graph signal processing sandryhaila2013discrete ; chen2015discrete ; duvenaud2015convolutional . On the spectral domain, the operator $Δ$ is a real-valued symmetric semidefinite matrix and the graph convolution is parameterized by a learnable filter $g_{θ}$

on the its eigenvalue matrix. Kipf et al.

kipf2016semi made assumptions of the largest eigenvalue (i.e.,

λ_{m a x} = 2

) and simplified it with two-order Chebyshev expansion,

g_{θ} ⋆ x \approx θ_{0} x + θ_{1} (Δ - I) x \approx θ (I + D^{- \frac{1}{2}} A D^{- \frac{1}{2}}) x .

(3)

A multi-layer graph convolutional network (GCN) is formulated as the following layer-wise propagation rule ( $σ (\cdot)$

is an activation function, e.g., ReLU):

X^{(l + 1)} = σ ({~ D}^{- \frac{1}{2}} (I + A) {~ D}^{- \frac{1}{2}} X^{(l)} W^{(l)})

(4)

where ${~ D}^{- \frac{1}{2}} (I + A) {~ D}^{- \frac{1}{2}} \leftarrow I + D^{- \frac{1}{} 2} A D^{- \frac{1}{2}}$ is the renormalization trick, $X^{(l)}$ and $W^{(l)}$ are the layer-wise feature and parameter matrices, respectively.

3 GCN as Two-Step Optimization

In this section, we re-interpret a GCN as a two-step optimization problem, where STEP1 is to minimize $L_{r e g}$ by viewing ${W^{(l)}}$ as constants and ${X^{(l)}}$ as parameters while STEP2 minimizes $L_{0}$ by updating ${W^{(l)}}$ . Overall, the GCN architecture is interpreted as a layer-wise combination of MLP architecture and the gradient descent algorithm of minimizing $L_{r e g}$ . In the meantime, the training process of parameters is entirely inherited from MLP, which aims at only minimizing $L_{0}$ . Let us first discuss a gradient descent algorithm for minimizing the form $\frac{1}{2} Tr (X^{⊤} Δ X)$ .

3.1 Gradient Descent for Trace Optimization.

Problem Definition.

Given the Laplacian operator $Δ \in R^{n \times n}$ , we consider to minimize the trace on feature domain $X \in R^{n \times m}$ , where $m$ is the input dimension. To prevent the trivial solution $X = 0 \in R^{n \times m}$ , we consider the energy constraint on $X$ , i.e., $∥ X ∥_{F}^{2} = c_{1} \in R^{+}$ . The trace optimization problem is:

min \frac{1}{2} Tr (X^{⊤} Δ X), s u b j e c t t o c o n s t . ∥ X ∥_{F}^{2},

(5)

where $∥ X ∥_{F}^{2}$ denotes the Forbenius-norm of $X$ . To solve this, We equivalently transform the optimization problem into the Reyleigh Quotient form $R (X)$ , which is,

min R (X) = \frac{\frac{1}{2} Tr (X^{⊤} Δ X)}{∥ X ∥_{F}^{2}} = \frac{\frac{1}{2} Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}, s u b j e c t t o c o n s t . ∥ X ∥_{F}^{2} .

(6)

It is obvious that $R (X)$ is scaling invariant on $X$ , i.e., $\forall c_{2} \neq 0 \in R$ , $R (X) = R (c_{2} \cdot X)$ .

One-step Improvement.

Suppose an initial guess is $X$ , one-step of trace optimization aims at finding a better guess $X^{''}$ , which satisfies $R (X^{''}) \leq R (X)$ and $∥ X^{''} ∥_{F}^{2} = c_{1}$ . Our strategy is first view the problem as unconstrainted optimization on $R (X)$ and update the guess $X$ to $X^{'}$ by gradient descent. Then we rescale $X^{'}$ and reach the improved guess $X^{''}$ , which meets the norm constraint.

Given the initial guess $X$ , we move against the derivative of $R (X)$ by the learning rate $η = Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)})$ and reach an intermediate solution $X^{'}$ in the unconstrainted space:

	$\nabla_{X}$	$= \frac{\partial R (X)}{\partial X} = \frac{1}{2} \frac{\partial \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}}{\partial X} = \frac{(Δ - I \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}) X}{Tr (X^{⊤} X)},$		(7)
	$X^{'}$	$= X - η \nabla_{X} = \frac{(2 - Δ) X}{2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}} = \frac{(I + D^{- \frac{1}{2}} A D^{- \frac{1}{2}})}{} 2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)} X .$		(8)

Immediately, we get $R (X^{'}) \leq R (X)$ (details in Appendix B). Then, we rescale $X^{'}$ (to achieve the improved guess $X^{''} = c_{2} \cdot X^{'}$ , which naturally satisfies $R (X^{''}) = R (X^{'})$ ) by a constant $c_{2} \in R^{+}$ to meet the norm constraint, such that $∥ X^{''} ∥_{F}^{2} = c_{1}$ . Therefore, the improved guess satisfies $R (X^{''}) = R (X^{'}) \leq R (X)$ and has the following form,

X^{''} = c_{2} \cdot X^{'} \propto (I + D^{- \frac{1}{2}} A D^{- \frac{1}{2}}) X .

(9)

Note that, if we conduct the trace optimization algorithm enough times, the optimal solution will finally be proportional to the largest eigenvector of

Δ

, which causes oversmoothing.

3.2 Layer-wise Propagation and Optimization

We introduce the trace optimization solution into the layer-wise propagation of MLP. Given the node set $V$ , features $X = {x_{1}, x_{2}, \dots, x_{n}}$ and a labelled set $T \subset V$ , a label mapping, $f_{{W^{(l)}}} : X \mapsto Y$ , is usually a deep neural network, which could be tailored according to the practical applications. For example, $f (\cdot)$

could be a convolutional neural network (CNN) for image recognition or a recurrent neural network (RNN) for language processing. In this scenario, we first consider a simple multi-layer perceptron (MLP). The forward propagation rule of an MLP is given by,

X^{(l + 1)} = σ (X^{(l)} W^{(l)}),

(10)

where $W^{(l)}$ and $X^{(l)}$ are layer-wise parameters and inputs.

STEP1: minimizing $L_{r e g}$ in Forward Propagation.

Let us fix parameters ${W^{(l)}}$ and consider $X^{(l)}$ , i.e., the output of $(l - 1)_{t h}$ layer, as an initial guess of the trace optimization problem. We know from Sec. 3.1 that through one-step gradient descent, Eqn. (9) will find an improved guess for minimizing $L_{r e g}$ :

X^{(l)^{''}} = c_{2} \cdot X^{(l)^{'}} = {X^{(l)} - η \nabla_{X} ∣ ∣}_{X = X^{(l)}}^{(l)} \propto (I + D^{- \frac{1}{2}} A D^{- \frac{1}{2}}) X^{(l)} .

(11)

We plug this new value into Eqn. (10) and immediately reach the same convolutional propagation rule, $f^{(l)} : X^{(l)} \mapsto X^{(l + 1)}$ , as Kipf et al. kipf2016semi (before applying the renormalization trick),

X^{(l + 1)} = σ ((I + D^{- \frac{1}{2}} A D^{- \frac{1}{2}}) X^{(l)} W^{(l)}),

(12)

where the constant scalar in Eqn. (11) is absorbed into parameter matrix $W^{(l)}$ . Therefore, a GCN forward propagation is essentially applying STEP1 layerwise in the forward propagation of an MLP, which is formulated as a composition of mappings $f = f^{(L - 1)} \circ \dots \circ f^{(0)}$ on initial feature $X^{(0)}$ .

STEP2: minimizing $L_{0}$ in Back Propagation.

After forward propagation, the cross-entropy loss $L_{0}$ over the labelled set is calculated. In this procedure, we then conversely view ${X^{(l)}}$ as constants and ${W^{(l)}}$ as parameters of the MLP and conduct standard back-propagation algorithm.

3.3 GCN: combining STEP1 and STEP2

In essence, STEP1 essentially defines a combined architecture, where the layer-wise propagation is adjusted by an additional step of trace optimization. In STEP2, under that architecture, the optimal ${W^{(l)}}$ is learned and a low-dimension $f (X)$ is reached with respect to $L_{0}$ explicitly and $L_{r e g}$ implicitly, after standard loss back-propagation. We express it as a two-step optimization,

STEP1: min {X^{(l)}} L_{r e g} ({X^{(l)}} ∣ {W^{(l)}})      w h y G C N a r c h i t e c t u r e ? and STEP2: min {W^{(l)}} L_{0} ({W^{(l)}})      t r a i n t h e a r c h i t e c t u r e . .

(13)

In this section, the learning rate $η = Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)})$ is specially chosen, and it satisfies $η \in (0, \infty)$ since $X^{⊤} X$ is semi-definite and $\frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}$ is smaller than the largest eigenvalue of $Δ$ , which is smaller than $2$ . In the experiment section, we reveal that $η$ is related to the weight of neighbor information averaging. We further test different $η$ and provide more insights on graph convolution operators. In the following sections, we use $A_{s y m}$ to denote the convolutional operator ${~ D}^{- \frac{1}{2}} (I + A) {~ D}^{- \frac{1}{2}}$ and use $A_{r w}$ for the random walk form ${~ D}^{- 1} (I + A)$ .

4 The Over-smoothing Issue

The recent successes in applying GNNs are largely limited to shallow architectures (e.g., 2-4 layers). Model performance decreases when adding more intermediate layers. Summarized in zhao2019pairnorm , there are three possible contributing factors: (i) overfitting due to increasing number of parameters (one matrix $W^{(l)}$ per layer); (ii) gradient vanishing/exploding; (iii) oversmoothing due to Laplacian smoothing. The first two points are common to all deep neural networks. The issue of oversmoothing is therefore the focus in this work. We show that deep GCNs can learn anti-oversmoothing by nature, but overfitting is the major cause of performance drop.

In this section, we first recall the commonly discussed oversmoothing problem (in SGC). Based on the re-formulation in Sec. 3.3, we then show how deep GCNs have the ability to learn anti-oversmoothing. Further, we propose an easy but effective mean-subtraction trick to speed up anti-oversmoothing, which accelerates the convergence in training deep GCNs.

4.1 Over-smoothing in SGC.

Oversmoothing means that node representations become similar and finally go indistinguishable after multi-layer graph convolution. Starting from the random walk theory, the analysis for oversmoothing is usually done on a connected, un-directed and non-bipartite graph. The issue is discussed in li2018deeper ; xu2018representation ; li2019deepgcns ; rong2019dropedge ; zhao2019pairnorm and mainly on the $L$ -layer linear SGC wu2019simplifying .

SGC was proposed in wu2019simplifying , with the hypothesis that the non-linear activation is not critical while the majority of the benefit arises from the local averaging $A$ . The authors directly remove the activation function and proposed a linear “ $L$ -layer” model,

f (X) = A (A (\dots) W^{(L - 1)}) W^{(L)} = A^{L} X W .

(14)

where $\prod W^{(l)}$ has collapsed into a single $W$ . This model explicitly disentangles the dependence of STEP1 and STEP2. We similarly formulate the SGC model in the form of two-step optimization,

STEP1: min X L_{r e g} (X) and STEP2: min W L_{0} (W)

(15)

Theorem 1.

Given any random signal $x \in R^{n}$ and a symmetric matrix $A \in R^{n \times n}$ , the following property $lim k \to \infty A^{k} x \propto u_{1}$ holds almost everywhere on $x$ , where $A$ has non-negative eigenvalues and $u_{1}$ is the eigenvector associated with the largest eigenvalue of $A$ .

For two widely used convolution operators $A_{s y m}$ and $A_{r w}$ , they have the same dominant eigenvalue $λ_{1} = 1$ with eigenvectors ${~ D}^{\frac{1}{2}} 1$ and 1, respectively. From the two-step optimization form, SGC is essentially conducting gradient descent algorithm $L$ times in STEP1. According to Theorem 1 (see proofs in Appendix D), if $L$ goes to infinity, then each output feature channel will become $lim k \to \infty A_{s y m}^{k} x \propto {~ D}^{\frac{1}{2}} 1 \in R^{n}$ or $lim k \to \infty A_{r w}^{k} x \propto 1 \in R^{n}$ , which means oversmoothing. In STEP2, SGC model will seek to minimize $L_{0}$ on the basis of the oversmoothed features. The independence between STEP1 and STEP2 accounts for the performance drop in deep SGC.

4.2 Anti-oversmoothing in GCN.

On the contrary, the result of $L_{r e g}$ in STEP1 is dependent on ${W^{(l)}}$ , i.e., also on STEP2, in GCN. In fact, after STEP1, node representation in GCN will be oversmoothed to some extent as well. However, during STEP2, GCN will learn to update layer-wise ${W^{(l)}}$ and make node features separable, such that $L_{0}$ will be minimized, during which, the effect of STEP1 (minimizing $L_{r e g}$ , i.e., making node features inseparable) will be mitigated and $L_{r e g}$ actually increases implicitly. In essence, the dependency of $L_{r e g}$ enables GCNs to do anti-ovesmoothing during STEP2.

We demonstrate on the Karate club dataset: this graph has 34 vertices of 4 classes (the same labeling as kipf2016semi ; perozzi2014deepwalk

) and 78 edges. A 32-layer GCN (deep enough for this demo dataset), with 16 hidden units in each layer, is considered. The model is trained on identity feature matrix with basic residual connection

he2016deep . The training set

T

consists of two labeled examples per class. After 1000 epochs, the model achieves

96.15 %

accuracy in the testing samples. We present the feature-wise smoothing by layer in Fig. 2.a and node-wise smoothing by layer in Fig. 2

.b. The y-axis score of the first two figures are calculated by cosine similarity. We also calculate and present

L_{0}

as well as

L_{r e g}

for each training epoch in Fig. 2.c. More details of the demo are given in Appendix F.

[width=5.2in]fig/karate.pdf

Figure 2: Feature-wise and Node-wise Smoothing during GCN Training on Karate

From the demonstrations, we observe that without training (blue curves in Fig. 2.a and Fig. 2.b and the beginning of Fig. 2.c), the issue of oversmoothing do exist for deep GCNs. Because forward propagation mixes up features layer-by-layer. However, this issue is automatically addressed during training and with more training epochs, feature-wise smoothing and node-wise smoothing are gradually mitigated. The effect is more obvious when referring to Fig. 2.c, where $L_{r e g}$ actually increases during the training (small $L_{r e g}$ indicates oversmoothing). The gradual increase in $L_{r e g}$ demonstrates that GCNs have the ability to learn to anti-oversmooth by nature. Then, the next question comes: what is the real cause of performance drop in deep GCNs? Our answer is overfitting. We give practical support in the experimental section.

4.3 Mean-subtraction: an Accelerator

Although deep GCNs could learn anti-oversmoothing naturally, another practical problem appears that the convergence of training deep GCNs is extremely slow (sometimes may not converge). This issue has not been explored extensively in the literature. In this work, we present a mean-subtraction trick to accelerate the training of deep GCNs, which theoretically magnifies the effect of Fiedler vector. PairNorm zhao2019pairnorm also includes a mean-subtraction step, however, their purpose is to simplify derivation. This section provides more insights and motivation to use mean-subtraction.

There are primarily two reasons to use mean-subtraction: (i) deep neural network classifiers are discriminators that draw a boundary between classes. Therefore, the mean feature of the entire dataset (a DC signal) does not help with the classification, whereas the components away from the center (the AC signal) matters; (ii) layer-wise

mean-subtraction will eliminate the dominant eigen component (

{~ D}^{\frac{1}{2}} 1

1

) and actually magnifies the Fiedler vector (the eigenvector associated with the second smallest eigenvalue of

Δ

), which reveals important community information and graph conductance chen2017supervised . This helps to set an initial graph partition and spends up model training (STEP2).

We start with one of the most popular convolution operator $A_{r w}$ and its largest eigenvector $u_{1} = 1 \in R^{n}$ . Given any non-zero $x \in R^{n}$ , the mean-subtraction gives,

x_{n e w} \leftarrow x - x_{m e a n} = x - \frac{1 1^{⊤} x}{n} = x - ⟨ x,_{1} ⟩ \cdot_{1}

(16)

where $_{1} = \frac{u_{1}}{∥ u_{1} ∥}$ . Eqn. (16) reveals that mean-subtraction reduces the components aligned with ${u_{1}}$ -space. This is exactly a step of numerical approximation of the Fiedler vector, which sets the initial graph partition (demonstration in Appendix F) and makes the feature separable. For $A_{s y s}$ , the formulation could be adjusted by a factor of ${~ D}^{\frac{1}{2}}$ (refer to Appendix E).

5 Experiments

In this section, we present experimental evidence to answer the following three questions: (i) whether oversmoothing is an issue in deep GCNs and why? (ii) How to stabilize and accelerate the training of deep GCNs? (iii) Does the learning rate $η$ matter? How about changing them? We also provide more insights and draw useful conclusions for the practical usage of GCN models and its variants.

The experiments show the performance of deep GCNs on the semi-supervised node classification tasks. All the deep models (with more than 3 hidden layers) are implemented with basic skip-connection kipf2016semi ; he2016deep . Since skip-connection (also called residual connection) are necessary in deep architectures, we do not consider them as new models. Three benchmark citation networks (Cora, Citeseer, Pubmed) are considered. We follow the same experimental settings from yang2016revisiting and show the basic statistics of datasets in Table. 1. All the experiments are conducted 20 times and mainly finished in a Linux server with 64GB memory, 32 CPUs and a single GTX-2080 GPU.

Dataset	Nodes	Edges	Features	Class	Label rate
Cora	2,708	5,429	1,433	7	0.052
Citeseer	3,327	4,732	4,732	6	0.036
Pubmed	19,717	44,338	500	3	0.003

Table 1: Overview of Citation Network Statistics

5.1 Overfitting in Deep GCNs

The performance of GCNs is known to decrease with increasing number of layers, for which, a common explanation is “oversmoothing”. In Sec. 4, we contradict this thesis and conjecture instead that overfitting is the major reason for the drop of performance in deep GCNs; we show that deep GCNs actually learn anti-oversmoothing. In this section, we provide evidence to support our conjecture.

Performance vs Depth.

We first evaluate vanilla GCN models (with residual connection) with $1 \sim 50$ hidden layers on Cora, Citeseer and Pubmed. The results of training and test accuracy are reported in Fig. 3.

Form Fig. 3, we know immediately that test accuracy drops in the beginning (1-4 layers) and then remains stable even as the model depth increases, which means the increasing number ( $\geq 4$ ) of hidden layers does not hurt model performance. Thus, oversmoothing is not the reason. From 2 to 3 or 3 to 4 layers, we notice that these is a big rise in training accuracy (up to $\geq 99 %$ ) and simultaneously a big drop in training loss (to $\leq 0.1$ ) and test accuracy consistently on the three datasets. This is more consistent with overfitting.

[width=4.9in]fig/threesets.pdf

Figure 3: Vanilla Deep GCNs on Three Citation Networks

Deep GCNs Learn Anti-oversmoothing.

We recall that in Sec. 4.2, we show from an optimization perspective that the dependency of $L_{r e g}$ on ${W^{(l)}}$ allows the network to learn anti-oversmoothing. To verify our theory, we compare SGC and GCN on Cora and Pubmed with various depth. To make it clear, SGC is actually a linear model, the depth $L$ means the number of graph convolution $S^{L}$ . Model performance in both training and test is shown in Fig. 4.

It is interesting that the accuracy of SGC decreases rapidly oono2019graph with more graph convolutions either for training or test. This is a strong indicator of oversmoothing, because node features converge to the same stationary due to the effect of STEP1 (specified in Theorem 1). The performance of the GCN model is not as good as SGC soon after 2 layers because of overfitting, but it stabilizes at a high accuracy even as the model goes very deep, which again verifies that GCNs naturally have the power of anti-oversmoothing.

[width=5.5in]fig/GCNvsSGC2.pdf

Figure 4: Comparison of Deep GCN and Deep SGC on Cora and Pubmed

5.2 Mean-Subtraction

To facilitate the training of deep GCNs, we proposed mean-subtraction in Sec. 4.3. In this section, we evaluate the efficacy of the mean-subtraction trick and compare it with vanilla GCNs kipf2016semi , PairNorm zhao2019pairnorm and the widely used BatchNorm ioffe2015batch

in the deep learning area. The four models have same settings, such as the number of layers (64), learning rate (0.01), and hidden units (16). Mean-subtraction is to subtract the mean feature value before each convolution layer (and PairNorm further re-scales the feature by variance). They do not include additional parameters. BatchNorm adds more parameters for each layer and learns to whiten the input of each layer. The experiment is conducted on

Cora with

500

training epochs.

[width=4.5in]fig/mean.pdf

Figure 5: Comparison of Effective Tricks to Train Deep GCNs (20 trials)

Fig. 5 reports the training and test curve for the four model variants. GCN and GCN+BatchNorm perform similarly, which means BatchNorm does not help substantially in training deep GCNs (at least for Cora). GCN+mean-subtraction and GCN+PairNorm give fast and stable training/test convergence. However, we could tell from the training curve that the PairNorm trick seems to suffer a lot from overfitting, leading to a drop in test accuracy. In sum, mean-subtraction not only speeds up the model convergence but also retains the same expressive power. It is an ideal trick for training deep GCNs.

5.3 Performance vs Learning Rate $η$

In Sec. 3, we choose the learning rate $η = Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)})$ . However, a different learning rate does lead to different weights $w (η)$ of neighbor information aggregation (we show that $w (η)$ is a monotonically increasing function in Appendix C). There are also some efforts on trying different ways to aggregate neighbor information velivckovic2017graph ; hamilton2017inductive ; rong2019dropedge ; chen2018fastgcn . In this section, we consider the form “ $I + w (η) A$ ” with $w (η) \in [0, \infty)$ and exploit a group of convolution operators by their normalized version. GCN with normalized $I + w (η) A$ is named as $η$ -GCN. We evaluate this operator group on Cora and list the experimental results in Table. 2

Accuracy (%)		$w (η)$ =0	0.1	0.2	0.5	1.0	2	5	10	20	50	100
2-layer	training	92.66	95.67	96.32	96.05	95.33	94.54	93.44	93.30	92.82	92.86	92.98
2-layer	test	52.66	71.62	75.64	78.26	79.72	79.06	78.34	78.14	77.93	77.67	77.58
32-layer	training	95.02	99.49	99.58	99.35	98.69	98.10	98.84	98.83	98.81	98.76	98.83
32-layer	test	39.93	72.53	73.59	73.65	74.03	75.11	74.16	75.08	75.49	74.64	74.74

Table 2: Performance vs Neighbor Averaging Weight (2-layer and 32-layer)

We conclude that when $w (η)$ is small (i.e., $η$ is small), which means the gradient of $L_{r e g}$ does not contribute much to the end effect, $η$ -GCN is more of a DNN. As $w (η)$ increases, a significant increase in model performance is initially observed. When $w (η)$ exceeds some threshold, the accuracy saturates, remaining high (or maybe decreases slightly in shallow models, i.e., $2$ -layer) even as we increase $w (η)$ substantially. We conclude that for the widely used shallow GCNs, the common choice of weight $w (η) = 1$ , which means a learning rate, $η = Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)})$ , is large enough to include the gradient descent effect and small enough to avoid the drop in accuracy. To find the best weight of neighbor averaging, further inspection is needed in future work.

6 Conclusion

We reformulate GCNs from an optimization perspective by plugging the gradient of graph regularizer into a standard MLP. From this formulation, we revisit the commonly discussed “oversmoothing issue” in deep GCNs and provide a new understanding: deep GCNs have the power to learn anti-oversmoothing by nature, but overfitting is the real cause of performance drop when the model goes deep. We further propose a cheap but effective mean-subtraction trick to accelerate the training of deep GCNs. Extensive experiments are presented to verify our theory and provide more practical insights.

References

(1) Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec. Hierarchical graph representation learning with differentiable pooling. In Advances in neural information processing systems, pages 4800–4810, 2018.
(2) Petar Veličković, William Fedus, William L Hamilton, Pietro Liò, Yoshua Bengio, and R Devon Hjelm. Deep graph infomax. arXiv preprint arXiv:1809.10341, 2018.
(3) Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint arXiv:1609.02907, 2016.
(4) Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst. Convolutional neural networks on graphs with fast localized spectral filtering. In Advances in neural information processing systems, pages 3844–3852, 2016.
(5) 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, 2019.
(6) Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. arXiv preprint arXiv:1710.10903, 2017.
(7) Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. In Advances in neural information processing systems, pages 1024–1034, 2017.
(8) Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl. Neural message passing for quantum chemistry. In
Proceedings of the 34th International Conference on Machine Learning-Volume 70
, pages 1263–1272. JMLR. org, 2017.
(9) Qi Liu, Maximilian Nickel, and Douwe Kiela. Hyperbolic graph neural networks. In Advances in Neural Information Processing Systems, pages 8228–8239, 2019.
(10) Johannes Klicpera, Aleksandar Bojchevski, and Stephan Günnemann. Predict then propagate: Graph neural networks meet personalized pagerank. arXiv preprint arXiv:1810.05997, 2018.
(11) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful are graph neural networks? arXiv preprint arXiv:1810.00826, 2018.
(12) Kenta Oono and Taiji Suzuki. Graph neural networks exponentially lose expressive power for node classification. arXiv preprint cs.LG/1905.10947, 2019.
(13) Andreas Loukas. What graph neural networks cannot learn: depth vs width. arXiv preprint arXiv:1907.03199, 2019.
(14) Nima Dehmamy, Albert-László Barabási, and Rose Yu. Understanding the representation power of graph neural networks in learning graph topology. In Advances in Neural Information Processing Systems, pages 15387–15397, 2019.
(15) Guohao Li, Matthias Muller, Ali Thabet, and Bernard Ghanem. Deepgcns: Can gcns go as deep as cnns? In
Proceedings of the IEEE International Conference on Computer Vision
, pages 9267–9276, 2019.
(16) Qimai Li, Zhichao Han, and Xiao-Ming Wu.
Deeper insights into graph convolutional networks for semi-supervised learning.
In
Thirty-Second AAAI Conference on Artificial Intelligence
, 2018.
(17) Keyulu Xu, Chengtao Li, Yonglong Tian, Tomohiro Sonobe, Ken-ichi Kawarabayashi, and Stefanie Jegelka. Representation learning on graphs with jumping knowledge networks. arXiv preprint arXiv:1806.03536, 2018.
(18) Yu Rong, Wenbing Huang, Tingyang Xu, and Junzhou Huang. Dropedge: Towards deep graph convolutional networks on node classification. In International Conference on Learning Representations, 2019.
(19) Lingxiao Zhao and Leman Akoglu. Pairnorm: Tackling oversmoothing in gnns. arXiv preprint arXiv:1909.12223, 2019.
(20) Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps for dimensionality reduction and data representation. Neural computation, 15(6):1373–1396, 2003.
(21) Zhilin Yang, William W Cohen, and Ruslan Salakhutdinov. Revisiting semi-supervised learning with graph embeddings. arXiv preprint arXiv:1603.08861, 2016.
(22) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In
Proceedings of the IEEE conference on computer vision and pattern recognition
, pages 770–778, 2016.
(23) Joshua B Tenenbaum, Vin De Silva, and John C Langford. A global geometric framework for nonlinear dimensionality reduction. science, 290(5500):2319–2323, 2000.
(24) Xiaojin Zhu and Zoubin Ghahramani. Learning from labeled and unlabeled data with label propagation. Technical report, Carnegie Mellon University, 2002.
(25) Jason Weston, Frédéric Ratle, Hossein Mobahi, and Ronan Collobert. Deep learning via semi-supervised embedding. In Neural networks: Tricks of the trade, pages 639–655. Springer, 2012.
(26) Mikhail Belkin and Partha Niyogi. Laplacian eigenmaps and spectral techniques for embedding and clustering. In Advances in neural information processing systems, pages 585–591, 2002.
(27) Yu Chen, Lingfei Wu, and Mohammed J Zaki. Deep iterative and adaptive learning for graph neural networks. arXiv preprint arXiv:1912.07832, 2019.
(28) Rie K Ando and Tong Zhang. Learning on graph with laplacian regularization. In Advances in neural information processing systems, pages 25–32, 2007.
(29) Alexander J Smola and Risi Kondor. Kernels and regularization on graphs. In Learning theory and kernel machines, pages 144–158. Springer, 2003.
(30) Uri Shaham, Kelly Stanton, Henry Li, Boaz Nadler, Ronen Basri, and Yuval Kluger. Spectralnet: Spectral clustering using deep neural networks. arXiv preprint arXiv:1801.01587, 2018.
(31) Ulrike Von Luxburg.
A tutorial on spectral clustering.
Statistics and computing, 17(4):395–416, 2007.
(32) Aliaksei Sandryhaila and José MF Moura. Discrete signal processing on graphs. IEEE transactions on signal processing, 61(7):1644–1656, 2013.
(33) Siheng Chen, Rohan Varma, Aliaksei Sandryhaila, and Jelena Kovačević. Discrete signal processing on graphs: Sampling theory. IEEE transactions on signal processing, 63(24):6510–6523, 2015.
(34) David K Duvenaud, Dougal Maclaurin, Jorge Iparraguirre, Rafael Bombarell, Timothy Hirzel, Alán Aspuru-Guzik, and Ryan P Adams. Convolutional networks on graphs for learning molecular fingerprints. In Advances in neural information processing systems, pages 2224–2232, 2015.
(35) Bryan Perozzi, Rami Al-Rfou, and Steven Skiena. Deepwalk: Online learning of social representations. In Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 701–710, 2014.
(36) Zhengdao Chen, Xiang Li, and Joan Bruna. Supervised community detection with line graph neural networks. arXiv preprint arXiv:1705.08415, 2017.
(37) Sergey Ioffe and Christian Szegedy. Batch normalization: Accelerating deep network training by reducing internal covariate shift. arXiv preprint arXiv:1502.03167, 2015.
(38) Jie Chen, Tengfei Ma, and Cao Xiao. Fastgcn: fast learning with graph convolutional networks via importance sampling. arXiv preprint arXiv:1801.10247, 2018.

Appendix A $L_{r e g}$ and Spectral Clustering

Graph Regularizer $L_{r e g}$ .

$L_{r e g}$ is commonly formulated by Dirichlet energy, $L_{r e g} = \frac{1}{2} Tr (H^{⊤} Δ H)$ , where $f (\cdot)$ is a mapping from the input feature $X$ to low-dimensional representation $H = f (X)$ . To minimize $L_{r e g}$ , this paper adds constraint on the magnitude of $H$ , i.e., $∥ H ∥_{F}^{2} = C \in R$ , which gives,

min H \frac{1}{2} Tr (H^{⊤} Δ H), s u b j e c t t o c o n s t . ∥ H ∥_{F}^{2} .

(17)

Spectral Clustering.

Given a graph with binary adjacency matrix $A$ , a partition of node set $V$ into $k$ set could be written as $P_{1}, P_{2}, \dots, P_{k}$ in graph theory. For normalized spectral clustering, the $k$ indicator vectors is written as $h_{i} = (h_{i}^{1}, \dots, h_{i}^{n})$ , where $h_{i}^{j}$ represents the affiliation of node $j$ in class set $P_{i}$ and $v o l (P_{i}) = \sum_{v_{j} \in P_{i}} d_{j}$ is the volume.

h_{i}^{j} = {\begin{matrix} \frac{1}{\sqrt{v o l (P_{i})}}, & v_{j} \in P_{i} 0, & o t h e r w i s e \end{matrix}

(18)

The $H = [h_{i}^{j}]_{i = 1.. k, j = 1.. n}$ is a matrix containing these $k$ indicator vectors as columns. For each row of $H$ , there is only one non-empty entry, implying $h_{i}^{⊤} h_{j} = 0, \forall i \neq j$ . Let us revisit the Normalized Cut of a graph for a partition $P_{1}, P_{2}, \dots, P_{k}$ .

$Ncut (P_{1}, P_{2}, \dots, P_{k})$	$= k \sum i = 1 \frac{c u t (P_{i},_{i})}{v o l (P_{i})}$
	$= \frac{1}{2} k \sum i = 1 \frac{\sum_{v_{j} \in P_{i}, v_{t} \notin P_{i}} A_{j t} + \sum_{v_{j} \notin P_{i}, v_{t} \in P_{i}} A_{j t}}{\sqrt{v o l (P_{i})} \sqrt{v o l (P_{i})}}$
	$= \frac{1}{2} k \sum i = 1 ⎛ ⎝ \sum v_{j} \in P_{i}, v_{t} \notin P_{i} {(\frac{1}{\sqrt{v o l (P_{i})}})}^{2} + \sum v_{j} \notin P_{i}, v_{t} \in P_{i} {(\frac{1}{\sqrt{v o l (P_{i})}})}^{2} ⎞ ⎠$
	$= \frac{1}{2} k \sum i = 1 \sum j, t (h_{i}^{j} - h_{i}^{t})^{2} = \frac{1}{2} k \sum i = 1 h_{i}^{⊤} L h_{i} = \frac{1}{2} Tr (H^{⊤} L H) .$	(19)

Also, $H$ satisfies $H^{⊤} D H = I$ . When the discreteness condition is relaxed and $H$ is substitute by $H = D^{- \frac{1}{2}} U$ , the normalized graph cut problem (normalized spectral clustering) is relaxed into,

min U \frac{1}{2} Tr (U^{⊤} Δ U), s u b j e c t t o U^{⊤} U = I .

(20)

This is a standard trace minimization problem which is solved by the matrix the eigen matrix of $Δ$ . Compared to Eqn. (17), Eqn. (20) has a stronger contraints, which outputs the optimal solution irrelevant to the inputs (feature matrix $X$ ). However, Eqn. (17) only add constraints on the magnitude of $H$ , which balances the trade-off and will give a solution induced by both the eigen matrix of $Δ$ and the original feature $X$ .

Appendix B Reyleigh Quotient

Reyleigh Quotient.

The Reyleigh Quotient of a vector $x \in R^{m}$ is the scalar,

R (x) = \frac{x^{⊤} Δ x}{x^{⊤} x},

(21)

which is invariant to the scaling of $x$ . For example, $\forall c_{1} \neq 0 \in R$ , we have $R (x) = R (c_{1} \cdot x)$ . When we view $R (x)$ as a function on $m$ -dim variable $x$ , it has stationary points $x_{i}$ , where $x_{i}$ is the eigenvector of $Δ$ . Let us assume $Δ x_{i} = λ_{i} x_{i}$ , then the stationary value at point $x_{i}$ will be exactly the eigenvalue $λ_{i}$ ,

R (x_{i}) = \frac{x_{i}^{⊤} Δ x_{i}}{x_{i}^{⊤} x_{i}} = \frac{x_{i}^{⊤} λ_{i} x_{i}}{x_{i}^{⊤} x_{i}} = λ_{i} .

(22)

When $x$ is not an eigenvector of $Δ$ , the partial derivatives of $R (x)$ with respect to the vector coordinate $x_{j}$ is calculated as,

	$\nabla_{x_{j}} R (x)$	$= \frac{\partial R (x)}{\partial x_{j}} = \frac{\frac{\partial}{\partial x_{j}} (x^{⊤} Δ x)}{x^{⊤} x} - \frac{(x^{⊤} Δ x) \frac{\partial}{\partial x_{j}} (x^{⊤} x)}{(x^{⊤} x)^{2}}$
		$= \frac{2 (Δ x)_{j}}{x^{⊤} x} - \frac{(x^{⊤} Δ x) 2 x_{j}}{(x^{⊤} x)^{2}} = \frac{2}{x^{⊤} x} (Δ x - R (x) x)_{j}$		(23)

Thus, the derivative of $R (x)$ with respect to $x$ is collected as,

\nabla R (x) = \frac{2}{x^{⊤} x} (Δ x - R (x) x) .

(24)

Minimizing R(x).

Suppose $Δ = I - D^{- \frac{1}{2}} A D^{- \frac{1}{2}}$ is the normalized Laplacian matrix. Let us first consider to minimize $R (x)$ without any constraints. Since $Δ$

is a symmetric real-valued matrix, it could be factorized by Singular Value Decomposition,

Δ = U Λ U^{⊤} = s \sum i = 1 u_{i} λ_{i} u_{i}^{⊤}

(25)

where $s$ is the rank of $Δ$ and $0 = λ_{1} \leq \dots \leq λ_{s} < 2$ are the eigen values. For any non-zero vector $x$ , it is decomposed w.r.t. the eigen space of $Δ$ ,

x = ϵ + s \sum i = 1 c_{i} \cdot u_{i}

(26)

where ${c_{i}}$ is the coordinates and $ϵ$ is a component tangent to the eigen space spanned by ${u_{i}}$ . Let us consider the component of $x$ within the eigen space and discuss $ϵ$ later. Therefore, the Reyleigh Quotient $R (x)$ can be calculated by,

R (x) = \frac{x^{⊤} Δ x}{x^{⊤} x} = \frac{(\sum_{i = 1}^{s} c_{i} \cdot u_{i}^{⊤}) (\sum_{i = 1}^{s} u_{i} λ_{i} u_{i}^{⊤}) (\sum_{i = 1}^{s} c_{i} \cdot u_{i})}{(\sum_{i = 1}^{s} c_{i} \cdot u_{i}^{⊤}) (\sum_{i = 1}^{s} c_{i} \cdot u_{i})} = \frac{\sum_{i = 1}^{s} c_{i}^{2} λ_{i}}{\sum_{i = 1}^{s} c_{i}^{2}}

(27)

Recall the partial derivative of $R (x)$ w.r.t. $x$ in Eqn. (24). Think about to minimize $R (x)$ by gradient descent and always consider the learning rate (the same as what we used in the main text. The factor $\frac{1}{2}$ is from that the $R (x)$ in appendix does not have the scalar $\frac{1}{2}$ ),

η = \frac{1}{2} Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}) = \frac{1}{2} \cdot \frac{x^{⊤} x}{2 - \frac{x^{⊤} Δ x}{x^{⊤} x}} = \frac{1}{2} \cdot \frac{x^{⊤} x}{2 - R (x)} .

(28)

The initial $x$ is regarded as an starting point, and the next point $x^{'}$ is given by gradient descent,

x^{'} = x - η \nabla R (x) = x - \frac{1}{2} \cdot \frac{x^{⊤} x}{2 - R (x)} \frac{2}{x^{⊤} x} (Δ x - R (x) x) = \frac{2 I - Δ}{2 - R (x)} x .

(29)

The new Reyleigh Quotient value is,

R (x^{'}) = \frac{x^{' ⊤} Δ x^{'}}{x^{' ⊤} x^{'}} = \frac{{(\frac{2 I - Δ}{2 - R (x)} x)}^{⊤} Δ (\frac{2 I - Δ}{2 - R (x)} x)}{{(\frac{2 I - Δ}{2 - R (x)} x)}^{⊤} (\frac{2 I - Δ}{2 - R (x)} x)} = \frac{x^{⊤} (2 I - Δ) Δ (2 I - Δ) x}{x^{⊤} (2 I - Δ) (2 I - Δ) x} .

(30)

The eigen properties of $2 I - Δ$ could be derived from $Δ$ , where they have the same eigenvector, and any eigenvalue $λ$ of $Δ$ will adjust to be an eigenvalue $2 - λ$ of $2 I - Δ$ . Therefore, we do further derivation,

R (x^{'}) = \frac{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2} λ_{i}}{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}} .

(31)

So far, to get the ideal effect, a final check is needed: whether the Reyleigh Quotient does decrease after the gradient descent.

$R (x^{'}) - R (x)$	$= \frac{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2} λ_{i}}{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}} - \frac{\sum_{i = 1}^{s} c_{i}^{2} λ_{i}}{\sum_{i = 1}^{s} c_{i}^{2}}$
	$= \frac{(\sum_{i = 1}^{s} c_{i}^{2}) (\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2} λ_{i}) - (\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}) (\sum_{i = 1}^{s} c_{i}^{2} λ_{i})}{(\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}) (\sum_{i = 1}^{s} c_{i}^{2})}$
	$= \frac{\sum_{i, j} c_{i}^{2} c_{j}^{2} (λ_{i} - λ_{j}) (λ_{j} - λ_{i}) (4 - λ_{i} - λ_{j})}{(\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}) (\sum_{i = 1}^{s} c_{i}^{2})}$
	$= - \frac{\sum_{i, j} c_{i}^{2} c_{j}^{2} (λ_{i} - λ_{j})^{2} (4 - λ_{i} - λ_{j})}{(\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2}) (\sum_{i = 1}^{s} c_{i}^{2})} < 0$	(32)

Also, we show the asymptotic property of $R (x)$ in gradient descent,

lim t \to \infty R (x^{(t)}) = lim t \to \infty \frac{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2 t} λ_{i}}{\sum_{i = 1}^{s} c_{i}^{2} (2 - λ_{i})^{2 t}} = \frac{c_{2}^{2} λ_{2}}{c_{1}^{2}} \cdot lim t \to \infty {(\frac{2 - λ_{2}}{2 - λ_{1}})}^{2 t} = 0^{+}

(33)

where $x^{(t)}$ is the $t$ -th new point given by gradient descent. So far, we finish the proof of well-definedness of gradient descent with the $η = Tr (X^{⊤} X) / (2 - \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)})$ .

Remark 1.

In fact, as stated above, $R (x)$ is invariant to the scaling of $x$ , so we could scale $x$ on its magnitude, i.e., making $∥ x ∥ = c_{2} \in R^{+}$ as a constraint during the gradient descent iteration, all the properties and results still hold.

Remark 2.

In the main text, instead of using a vector $x$ , we use a feature matrix $X$ and define our Reyleigh Quotient by $R (X) = \frac{Tr (X^{⊤} Δ X)}{Tr (X^{⊤} X)}$ . In fact, different feature channels of $X$ could be viewed as independent vector signal $x_{i} \in R^{m}$ and for each channel, the same gradient descent analysis is applied. Therefore, we finish the detailed proof for our formulation in the main text, which is of the following form,

min R (X), s u b j e c t t o c o n s t . ∥ X ∥_{F}^{2} .

(34)

Appendix C Learning Rate $η$ and Neighbor Averaging Weight $w (η)$

We show the relation of learning rate $η$ and neighbor averaging weight $w (η)$ in this section (the derivation is in terms of the main text, so $\nabla R (x)$ does not have factor $2$ ).

x^{'} = x - η \nabla R (x) = x - η \frac{1}{x^{⊤} x} (Δ x - R (x) x) = (\frac{R (x) - 1}{x^{⊤} x} \cdot η + 1) (x + \frac{η \cdot D^{- \frac{1}{2}} A D^{- \frac{1}{} 2}}{(R (x) - 1) \cdot η + x^{⊤} x} x)

(35)

Thus, we have,

w (η) = \frac{η}{(R (x) - 1) \cdot η + x^{⊤} x} = \frac{1}{R (x) - 1 + \frac{x^{⊤} x}{η}},

(36)

According to the formulation, $w (η)$ is a monotonically increasing function on variable $η$ and is valid when $w (η) > 0$ . Therefore, when $R (x) \geq 1$ , the domain $η \in [0, \infty)$ and when $R (x) < 1$ (we know from Eqn. (33) that $R (x) \to 0^{+}$ ), the domain of the function is bounded, $η \in [0, \frac{x^{⊤} x}{1 - R (x)})$ .

Remark 3.

The choice of $η = Tr (X^{⊤} X) / (2 -$