Year: 2021

Abstract:
Graph convolution operator of the GCN model is originally motivated from a localized first-order approximation of spectral graph convolutions. This work stands on a different view; establishing a mathematical connection between graph convolution and graph-regularized PCA (GPCA). Based on this connection, the GCN architecture, shaped by stacking graph convolution layers, shares a close relationship with stacking GPCA. We empirically demonstrate that the unsupervised embeddings by GPCA paired with a 1or 2-layer MLP achieves similar or even better performance than many sophisticated baselines on semi-supervised node classification tasks across five datasets including Open Graph Benchmark. This suggests that the prowess of graph convolution is driven by graph based regularization. In addition, we extend GPCA to the (semi-)supervised setting and show that it is equivalent to GPCA on a graph extended with “ghost” edges between nodes of the same label. Finally, we capitalize on the discovered relationship to design an effective initialization strategy based on stacking GPCA, enabling GCN to converge faster and achieve robust performance at large number of layers.

heading: 1 INTRODUCTION
text: Graph neural networks (GNNs) are neural networks designed for the graph domain. Since the breakthrough of GCN (Kipf & Welling, 2017), which notably improved performance on the semisupervised node classification problem, many GNN variants have been proposed; including GAT (Veličković et al., 2018), GraphSAGE (Hamilton et al., 2017), DGI (Veličković et al., 2019), GIN (Xu et al., 2019), PPNP and APPNP (Klicpera et al., 2019), to name a few.
Despite the empirical successes of GNNs in both node-level and graph-level tasks, they remain not well understood due to limited systematic and theoretical analysis of GNNs. For example, researchers have found that GNNs, unlike their non-graph counterparts, suffer from performance degradation with increasing depth, their expressive power decaying exponentially in number of layers (Oono & Suzuki, 2020). Such behavior is only partially explained by the oversmoothing phenomenon (Li et al., 2018; Zhao & Akoglu, 2020). Another surprising observation shows that a Simplified Graph Convolution model, named SGC (Wu et al., 2019), can achieve similar performance to various more complex GNNs on a variety of node classification tasks. Moreover, a simple baseline that does not utilize the graph structure altogether performs similar to state-of-the-art GNNs on graph classification tasks (Errica et al., 2020). These observations call attention to studies for a better understanding of GNNs (NT & Maehara, 2019; Morris et al., 2019; Xu et al., 2019; Oono & Suzuki, 2020; Loukas, 2020; Srinivasan & Ribeiro, 2020). (See Sec. 2 for more on understanding GNNs.)
Toward a systematic analysis and better understanding of GNNs, we establish a connection between the graph convolution operator of GCN (and PPNP) and Graph-regularized PCA (GPCA) (Zhang & Zhao, 2012), and show the similarity between GCN and stacking GPCA. This connection provides a deeper understanding of GCN’s power and limitation. Empirically, we also find that GPCA performance matches that of many GNN baselines on benchmark semi-supervised node classification tasks. We argue that the simple GPCA should be a strong baseline in future. What is more, the unsupervised stacking GPCA can be viewed as “unsupervised GCN” and provides a straightforward, yet systematic way to initialize GCN training. We summarize our contributions as follows:
• Connection between Graph Convolution and GPCA: We establish the connection between the graph convolution operator of GCN (also PPNP) and the closed-form solution of graph-regularized PCA (GPCA) formulation. We demonstrate that a simple graph-regularized PCA paired with 1- or 2-layer MLP can achieve similar or even better results than state-of-the-art GNN baselines over
several benchmark datasets. We further extend GPCA to (semi-)supervised setting which can generate embeddings using information of labels, which yields better performance on 3 out of 5 datasets. The outstanding performance of simple GPCA supports that the prowess of GCN on node classification task comes from graph based regularization. This motivates the study and design of other graph regularization techniques in the future.
• GPCANET: New Stacking GPCA model: Capitalizing on the connection between GPCA and graph convolution, we design a new GNN model called GPCANET shaped by (1) stacking multiple GPCA layers and nonlinear transformations, and (2) fine-tuning end-to-end via supervised training. GPCANET is a generalized GCN model with adjustable hyperparameters that control the strength of graph regularization of each layer. We show that with stronger regularization, we can train GPCANET with fewer (1–3) layers and achieve comparable performance to much deeper GCNs.
• First initialization strategy for GNNs: Capitalizing on the connection between GCN and GPCANET, we design a new strategy to initialize GCN training based on stacking GPCA, outperforming the popular Xaiver initialization (Glorot & Bengio, 2010). We show that the GPCANET-initialization is extremely effective for training deeper GCNs, that significantly improves the convergence speed, performance, and robustness. Notably, GPCANET-initialization is general-purpose and also applies to other GNNs. To our knowledge, it is the first initialization method specifically designed for GNNs.
We open-source code at http://bit.ly/GPCANet. All datasets are public-domain.

heading: 2 RELATED WORK
text: Understanding GNNs. Our work concerns learning on a single graph, hence we limit discussion of related work to node-level GNNs. GCN’s graph convolution is originally motivated from the approximation of graph filters in graph signal processing (Kipf & Welling, 2017). NT & Maehara (2019) show that graph convolution only performs low-pass filtering on original feature vectors, and also state a connection between graph filtering and Laplacian regularized least squares. Motivated by the oversmoothing phenomenon of graph convolution, Oono & Suzuki (2020) theoretically prove that GCN can only preserve information of node degrees and connected components when the number of layers goes to infinity, under some conditions of GCN weights. Recently several papers revisited the connection of graph convolution to graph-regularized optimization problem (Li et al., 2019; Ma et al., 2020; Pan et al., 2021; Zhao & Akoglu, 2020; Zhu et al., 2021), which is originally discussed in graph signal processing (Shuman et al., 2013). More specifically, both Ma et al. (2020) and Zhu et al. (2021) relate graph-regularization optimization to several GNNs such as GCN (Kipf & Welling, 2017), APPNP (Klicpera et al., 2019), and GAT (Veličković et al., 2018). However, all previous work study these connections while ignoring the learnable parameters, which are essential for high-performance deep learning. Our work differs from these by establishing a stronger and closer connection to graph-regularized PCA that also takes learnable parameters into account.
Graph-regularized PCA. PCA and its variants are standard linear dimensionality reduction approaches. Several work extend PCA to graph-structured data, such as Graph-Laplacian PCA (Jiang et al., 2013) and Manifold-regularized Matrix Factorization (Zhang & Zhao, 2012). For other variants, see Shahid et al. (2016).
Stacking Models and Deep Learning. The connection between CNN and stacking PCA has been explored in PCANet (Chan et al., 2015), which demonstrated that the (unsupervised) simple stacking PCA works as well as supervised CNN over a large variety of vision tasks. The original PCANet is shallow and does not have nonlinear transformations, while PCANet+ (Low et al., 2017) overcomes these limitations and pushes the architecture much deeper. The idea of layerwise stacking for feature extraction is not new and was empirically observed to exhibit better representation ability in terms of classification. For a comprehensive review, we refer to Bengio et al. (2013).
Initialization. Traditionally, neural networks (NNs) were initialized with random weights generated from Gaussian distribution with zero mean and a small standard deviation (Krizhevsky et al., 2012). As training deeper NNs became extremely difficult due to vanishing gradient and activation functions, Glorot & Bengio (2010) provided a specific weight initialization formula, named Xavier initialization, based on variance analysis without considering activation function. Xavier initialization is widely used for any type of NN even today, and it is the main initialization strategy used for GNNs. Later, He et al. (2015) adapted Xavier initialization to ReLU activation by considering a multiplier. Taking another
direction, Saxe et al. (2013) analyzed the dynamics of training deep NNs and proposed random orthonormal initialization. Mishkin & Matas (2015) further improved orthonormal initialization for batch normalization (Ioffe & Szegedy, 2015). Different from these data-independent approaches, others (Krähenbühl et al., 2016; Seuret et al., 2017; Wagner et al., 2013) have employed datadependent techniques, like PCA, to initialize deep NNs. Although initialization has been widely studied for general NNs, no specific initialization has been proposed for GNNs. In this work, we propose a data-driven initialization technique (based on GPCA), specific to GNNs for the first time.

heading: 3 GRAPH CONVOLUTION AND GPCA

heading: 3.1 GRAPH CONVOLUTION
text: Consider a node-attributed input graph G = (V,E,X) with |V | = n nodes and |E| = m edges, where X ∈ Rn×d denotes the node feature matrix with d features. Broadly, graph convolution operation convolves the features (or representations) over the graph structure.
GCN. Similar to other neural networks stacked with repeated layers, GCN contains multiple graph convolution layers each of which is followed by a nonlinear activation. Let H(l) be the l-th hidden layer representation, then, each GCN layer performs
H(l+1) = σ(ÃsymH (l)W (l)) (1)
where Ãsym = D̃− 1 2 (A+ I)D̃− 1 2 denotes the n×n symmetrically normalized adjacency matrix with self-loops, D̃ is the diagonal degree matrix where D̃ii = 1 + ∑n j=1 Aij , W (l) depicts the l-th layer parameters (to be learned), and σ is the nonlinear activation function. Formally, graph convolution is parameterized with W and maps an input X to a new representation Z as
Z = ÃsymXW . (2)
PPNP. For PPNP (Klicpera et al., 2019), the features are first transformed by an MLP before convolving over the graph. Formally, the operation is revised as
Z = µ ( I − (1− µ)Ãsym )−1 MLPW (X) = ( I + αL̃ )−1 MLPW (X) (3)
where we replace µ with α = (1− µ)/µ, L̃ := I − Ãsym denotes the normalized graph Laplacian, and W depicts the learnable MLP parameters. As matrix inverse is expensive, an approximate version called APPNP that employs the power method (Golub & Van Loan, 1989) is often used in practice.

heading: 3.2 GRAPH-REGULARIZED PCA (GPCA)
text: As stated by Bengio et al. (2013), “Although depth is an important part of the story, many other priors are interesting and can be conveniently captured when the problem is cast as one of learning a representation.” GPCA is one such representation learning technique with a graph-based prior.
Standard PCA learns k-dimensional projections Z ∈ Rn×k of feature matrix X ∈ Rn×d, aiming to minimize the reconstruction error ∥X − ZWT ∥2F , (4) subject to W ∈ Rd×k being an orthonormal basis. GPCA extends this formalism to graph-structured data by additionally assuming either smoothing bases (Jiang et al., 2013) or smoothing projections (Zhang & Zhao, 2012) over the graph. In this work we consider the latter case where low-dimensional projections are smooth over the input graph G, where L̃ = I− Ãsym denotes its normalized Laplacian matrix. The objective formulation of GPCA is then given as
min Z,W
∥X − ZWT ∥2F + αTr(ZT L̃Z) s.t. WTW = I (5)
where α is a hyperparameter that balances reconstruction error and the variation of the projections over the graph. Note that the first part of Eq. equation 5, along with the constraint, corresponds to the objective of the original PCA, while the second part is a graph regularization term that aims to “smooth” the learned representations Z over the graph structure. As such, GPCA becomes the standard PCA when α = 0.
Similar to PCA, the problem (5) is non-convex but has a closed-form solution (Zhang & Zhao, 2012). Surprisingly, as we show, it has a close connection with the graph convolution formulation in Eq. equation 2. In the following, we give the GPCA solution and then detail its connection to graph convolution.
Theorem 3.1. GPCA with formulation shown in (5) has the optimal solution (Z∗,W ∗) following Z∗ = (I + αL̃)−1XW ∗ , and W ∗ = (w1,w2, ...,wk) (6)
where w1, ...,wk are the eigenvectors of XT (I+αL̃)−1X corresponding to the largest k eigenvalues.
Proof. The proof can be found in Appendix. A.1.

heading: 3.3 CONNECTION BETWEEN GCN AND GPCA
text: Let Φα := I + αL̃. The normalized Laplacian matrix L̃ has absolute eigenvalues bounded by 1, thus, all its positive powers have bounded operator norm. When α ≤ 1, Φ−1α can be decomposed into Taylor series as (I + αL̃)−1 = I − αL̃+ . . .+ (−α)tL̃t + . . .. The first-order truncated form (i.e. approximation) of the series is (I + αL̃)−1 ≈ I − αL̃ = (1− α)I + αÃsym . (7) When α = 1, the first-order approximation of Z∗ in Theorem 3.1 follows Z∗ ≈ ÃsymXW ∗ . (8) The (approximate) solution to GPCA in Eq. equation 8 matches the form of graph convolution operation in Eq. equation 2, with W ∗ plugged in as the eigenvectors of the matrix XTΦ−1α X . In other words, there exists some parameter W ∗ with which GCN becomes the first-order approximation of GPCA.
To reiterate, a key contribution of this work is to show that the graph convolution operation in GCN can be viewed as the first-order approximation of GPCA with α = 1 with a learnable W . Put differently, the first-order approximation of (unsupervised) GPCA with α = 1 can be viewed as a graph convolution with a fixed, data-driven W . In other words, Notably, for α < 1, Eq. equation 7 shows the connection between GPCA and graph convolution equipped with 1-step (scaled) residual connection.

heading: 3.4 CONNECTION BETWEEN PPNP AND GPCA
text: Replacing the MLP in Eq. equation 3 with a single linear layer without activation results in Z =( I + αL̃ )−1 XW , which has exactly the same formulation as the solution Z∗ in Theorem 3.1 Eq.
equation 6. The connection states that the graph convolution in PPNP can be viewed as the GPCA solution with a learnable W . Interestingly, the empirical performance improvement of PPNP over GCN (see Table 2 in Klicpera et al. (2019)) may be explained through these connections that they have to GPCA; where PPNP relates to the exact solution of GPCA while GCN is related to its (first-order) approximation.

heading: 3.5 SUPERVISED GPCA
text: The standard GPCA problem in (5) is unsupervised. Motiviated from LDA (Balakrishnama & Ganapathiraju, 1998) and PLS (Geladi & Kowalski, 1986), in this section we show how to extend it to the supervised setting, by learning embeddings that not only (1) provide good reconstruction and (2) vary smoothly over the graph structure, but also (3) highly correlate with the response variable(s). For simplicity of presentation, let z ∈ Rd be a 1-d embedding and Y denote the response matrix (in the general case of multiple responses). We write the additional, i.e. (3)rd objective above, as1
max z
[ corr(Y, z) ]T [ corr(Y, z) ] var(z) ≡ max
z zTY Y T z (9)
The form of equation 9 (See Appendix. A.3) and the variance-maximizing term var(z) are for mathematical convenience. Despite agnostic to labels, including var(z) is intuitive since an implicit objective of data projection (embedding) is to ensure that inherent variation in data is captured as much as possible. In general, we would aim to maximize the trace of ZTY Y TZ for multi-dimensional embeddings.
Interpretation. For semi-supervised node classification with c classes, let L ⊂ V denote the set of labeled nodes. For this task, Y ∈ {0, 1}n×c would encode the node labels where the v-th row of Y , denoted Yv, depicts the one-hot encoded label for each v ∈ L. For u ∈ V \L with unknown labels, Yu = 0, set as the c-dimensional all-zero vector. Then, (Y Y T )ij is simply equal to 1 when
1For the optimization to be well-posed, constraints on z are required, omitted for simplicity of presentation.
nodes i and j share the same label, and otherwise 0 (either because they have different labels or labels are unknown). This term simply enforces the representations Zi and Zj of two same-labeled nodes to be similar. In a sense, Y Y T adds “ghost” edges between the same-label nodes, further guiding the smoothness of their representations over this extended graph structure. We remark that earlier work (Gallagher et al., 2008) has heuristically introduced edges between same-label nodes to enhance a given graph for the node classification task. In this work, we have derived the theoretical underpinning for this strategy.
Supervised formulation. We have shown that requiring the embeddings to correlate with the known labels can be interpreted as additional smoothing over “ghost” edges between the same-label nodes in the graph. As such, we extend the GPCA problem in (5) to the (semi-)supervised setting as
min Z,W
∥X − ZWT ∥2F + αTr(ZT L̃sprZ) s.t. WTW = I ; (10)
where L̃spr = I − Ãspr , Ãspr = (1− β)Ãsym + βD− 1 2 (Y Y T )D− 1 2 (11)
In Eq. equation 11, β is an additional hyperparameter for trading-off the graph-based regularization (i.e. smoothing) due to the actual input graph edges versus the ones introduced through Y Y T between the nodes of the same label, and D is the diagonal matrix with Dii = ∑n j=1(Y Y T )ij .
Theorem 3.2. Supervised GPCA, as shown in (10) has the optimal solution (Z∗,W ∗) following Z∗ = (I + αL̃spr) −1XW ∗ , and W ∗ = (w1,w2, ...,wk) (12)
where w1, . . . ,wk are the top eigenvectors of the matrix XT (I +αL̃spr)−1X , equivalently XT ( (1+
α)I − [ α(1− β)Ãsym + αβD− 1 2Y Y TD− 1 2 ])−1 X , corresponding to the largest k eigenvalues.
Proof. The proof is similar to that of Theorem 3.1.

heading: 3.6 APPROXIMATION AND COMPLEXITY ANALYSIS
text: According to formulations in Theorems 3.1 and 3.2, obtaining Z∗∈Rn×k and W ∗∈Rd×k requires two demanding computations (1) the inverse of Φα = (I + αL) ∈ Rn×n, or in the supervised case Φα = (I+αL̃spr); and (2) top k eigenvectors of the matrix XTΦ−1α X ∈ Rd×d. Eigen-decomposition takes O(d3) (Pan & Chen, 1999), which is scalable as d is usually small. Computing matrix inverse, on the other hand, can take O(n3) and require O(n2) memory, which would be infeasible for very large graphs.
To reduce computation and memory complexity, we instead approximately compute F := ϕ−1α X , which is a common term for both Z∗ and W ∗. We can equivalently write
(I + αL)F = X =⇒ F + αF = αPF +X =⇒ F = α 1 + α PF + 1 1 + α X
for P = Ãsym in the unsupervised case and P = (1−β)Ãsym +βD− 1 2 (Y Y T )D− 1 2 when supervised.
Then, we can iteratively (with total T iterations) use the power method (Golub & Van Loan, 1989) to compute F as
F (t+1) ← α 1 + α PF (t) + 1 1 + α X (13)
where t ∈ {0, ..., T} depicts the iteration and F (0) ∈ Rn×d is initialized as X (or randomly). For the supervised case, PF (t) is computed through a series of (from right to left) matrix-matrix products. This avoids the explicit construction of matrix Y Y T in memory. Overall, solving for F takes O(T (m+ n)d) where m is the number of edges in the graph. The supervised case has an additional term O(Td|L|c) with c being the number of classes and |L| ≤ n be the number of labeled nodes, which can also be upper-bounded by O(T (m+ n)d) when treating c as constant.
Having solved for F , we perform the matrix-matrix product Z∗ = FW ∗ in O(ndk) and then the eigen-decomposition of XTF in O(d3 + nd2) = O(nd2) (for n ≥ d). Assuming O(d) = O(k), overall complexity for computing the 1-layer GPCA is given as O(Tmd+ Tnd + nd2), which is linear in the number of nodes and edges. Note that empirically we found 5 ≤ T ≤ 10 to be sufficient.
Algorithm 1 GPCANET Forward Pass and Pre-training 1: Input: graph G = (V,E,X), GPCA hyper-parameter(s) α (and β if supervised, β = 0 otherwise), #layers
L, hidden layer sizes {d1, . . . , dL}, activation function σ(·), #approximation steps T 2: Output: pre-set layer-wise parameters {W (1), . . . ,W (L)} 3: Initialize H(0) := X 4: for l = 1 to L do 5: Center H(l−1) by subtracting mean of row vectors 6: F ←− H(l−1) 7: for t = 1 to T do 8: PF ←− (1− β)ÃsymF + βD− 1 2 (Y Y T )D− 1 2F 9: F ←− α 1+α PF + 1 1+α H(l−1)
10: end for 11: W (l) ←− top dl eigenvectors of H(l−1) T F 12: H(l) ←− σ(FW (l)) 13: end for

heading: 4 GPCANET: A STACKING GPCA MODEL

heading: 4.1 GPCANET
text: Thus far, we drew a connection between the geometrically motivated, manifold-based GPCA and the graph convolution operation of deep NN based GCN. Next we leverage this connection to design a new model called GPCANET that takes advantage of the relative strengths of each paradigm; namely, GPCA’s ability to capture data variation and structure, and GCN’s ability to capture multiple levels of abstraction (i.e. high-level concepts) through stacked layers and non-linearity.
In a nustshell, GPCANET is a stacking of multiple (unsupervised or supervised) GPCA layers and nonlinear transformations, which shares the same architecture as a multi-layer GCN. It consists of two main stages: (1) Pre-training, which initializes the layer-wise parameters through closed-form GPCA solutions, and (2) End-to-end-training, which refines these parameters through end-to-end gradient-based minimization of a global supervised loss criterion at the output layer.
We remark that GPCANET is not the same as GCN, as each layer uses the formulation in Thm.s 3.1 and 3.2 (with approximation shown in Sec. 3.6). In fact, when α = 1 and β = 0, GPCANET is the GCN model initialized with GPCANET-initialization, which we discuss more in Sec. 4.2. In other words, GPCANET is a generalized GCN model with additional hyperparameters, α and β, controlling the strength of graph regularization based on the existing or “ghost” edges, respectively.
Forward Pass and Pre-training stage. During pre-training, weights of the l-th layer, denoted as W (l) ∈ Rdl−1×dl , are pre-set (i.e. initialized) as the leading dl eigenvectors of the matrix H(l−1) T Φ−1α H
(l−1),2 where H(l−1) is the representation as output by the (l − 1)-th layer (with H(0) := X), and Φα can be the unsupervised (I+αL) or the supervised (I+αL̃spr). The pre-training stage takes a single forward pass. Algo. 1 shows both forward pass during end-to-end-training and the pre-training procedure, where line 11 in blue is a step used only for pre-training.
Additional treatment for ReLU: Nonlinear transformations like ReLU improves model capacity, however at pre-training stage, it causes information loss as all negative values are truncated to 0. This hinders the advantage of using the leading dl eigenvectors to initialize the weights so as to convey maximum variance (i.e. information) to the next layers. To address this issue, we instead use the leading dl/2 eigenvectors {wi}dl/2i=1 and their negatives {−wi} dl/2 i=1 to initialize W
(l). Empirically we observe this always improves performance when using ReLU activation.
End-to-end training stage. Pre-training can be seen as an information-preserving initialization, as compared to an uninformative random initialization, after which we refine the layer-wise parameters via gradient-based optimization w.r.t. a supervised loss criterion at the output layer. Specifically for semi-supervised node classification, we perform an end-to-end training w.r.t. the cross-entropy loss on the labeled nodes. All parameters are updated jointly through backpropagation during this stage, with forward computation shown in Algo.1 (excluding line 11).
2If d(l) is greater than the number of eigenvectors, all eigenvectors are used, with additional vectors generated from random projection of eigenvectors.

heading: 4.2 GPCANET-INITIALIZATION FOR GCN
text: When we set α = 1, β = 0, and approximate the matrix inverse (I + αL)−1 via first-order truncated Taylor expansion as shown in Eq. equation 7 , GPCANET has the same architecture with GCN. As such, we can use the pre-training stage of GPCANET to initialize GCN with only minor modification. Specifically, we replace lines 6 through 10 in Algo. 1 with the following single line:
F ←− ÃsymH(l−1) (14) The modified initialization is for GCN and is driven by the mathematical connection between GPCANET and GCN that we established. We expect that adapting it for other GNNs is also possible although we do not pursue this direction here.

heading: 5 EXPERIMENTS
text: In this section we design extensive experiments to answer the following questions. (Q1) How does the simple, unsupervised and shallow GPCA compare to its multi-layer extension GPCANET, as well as to existing GNNs? (Q2) How does our extended, semi-supervised GPCA compare to the original, unsupervised GPCA? (Q3) Does GPCANET-initialization improve GCN accuracy and robustness?

heading: 5.1 EXPERIMENTAL SETUP
text: Datasets. We focus on semi-supervised node classification (SSNC) and use 5 benchmark datasets: First three, CORA, CITESEER, PUBMED (Sen et al., 2008), are relatively small (2K to 10K nodes) but widely-used citation graphs. For these we use the data splits in Kipf & Welling (2017). The others, ARXIV and PRODUCTS, are newest and much larger (100K to 2000K) node classification benchmarks from Open Graph Benchmark (Hu et al., 2020), for which we use the official data splits. Data statistics can be found in Appendix. A.4. Baselines. We compare (unsupervised & semi-supervised) GPCA and GPCANET to state-of-the-art (SOTA) GNNs, including GCN (Kipf & Welling, 2017), APPNP (Klicpera et al., 2019), GAT (Veličković et al., 2018), and GraphSAGE (G-SAGE) (Hamilton et al., 2017). Model configuration and training. For each dataset, we define a separate pool of values for the hyperparameters (HPs): learning rate, weight decay, number of layers, hidden size, dropout rate, and regularization trade-off terms α, β. For fair comparison, all models share the same HP pools during training. See Appendix. A.5 for HP configurations and other details.

heading: 5.2 Q1: PERFORMANCE OF (UNSUPERVISED) GPCA AND GPCANET
text: GPCA. Having proved the mathematical connection between GPCA, GCN, and PPNP, we expect unsupervised GPCA (β = 0) to generate comparable representations. We perform GPCA with different α ∈ {1, 5, 10, 20, 50} (Appendix. Table 5) to obtain node representations and pass those to a 1- or 2-layer MLP. We compare to GCN, APPNP, as well as other GNNs; GAT and G-SAGE.
The performance results are given in Table 1. Due to the scale of the largest two datasets, ARXIV and PRODUCTS, we list the reported performance at OGB-leaderboard3 (depicted by ∗) for G-SAGE on both datasets, and that of (Cluster-)GAT on PRODUCTS.
We find that the simple 1-layer GPCA paired with MLP performs consistently better than the multilayer GCN model across all 5 datasets. GPCA’s performance is also comparable to or better than other SOTA GNNs. This is quite notable, given that GPCA is not only shallow but also unsupervised, whereas all other baselines are trained end-to-end, and with the exception of APPNP, they exhibit a multi-layer architecture. By carefully looking at the performance of GPCA with varying α (see Appendix. A.6), we find that different datasets have different best selected α∗ (in Table 1 top to bottom: α∗ = {50, 5, 10, 20, 20}) but in general a relatively larger α (compared to graph convolution of GCN that is equivalent to α = 1) is preferable for all datasets. Larger α implies stronger graphregularization on the representations. The outstanding performance of the simple GPCA empirically confirms that the power of GNNs on the SSNC problem is mainly driven by graph regularization.
GPCANET. Compared to the 1-layer GPCA, GPCANET has a deeper architecture along with nonlinear activation function. Moreover, it employs hyperparameter α at every layer to control the degree of graph regularization. As each graph convolution has fixed level of graph regularization, one may hypothesize that increasing the number of layers (L) corresponds to increasing the degree
3https://ogb.stanford.edu/docs/leader_nodeprop/
of graph regularization. We empirically test this hypothesis using GPCANET, by varying both L (2 to 10) and α (0.1 to 10) to show their connection (hidden size is fixed as 128). The result is shown in Figure 1. The diagonal pattern (in dark blue) empirically suggests that increasing the number of layers has the same effect as increasing graph regularization via α.
The corresponding interaction between α and number of layers suggests that we can train a GPCANET with fewer number of layers yet achieve similar regularization by increasing α. Such a shallow model that in fact behaves like a deep one has the advantage of less memory requirement and faster training due to fewer parameters.
To this end, we train 1–3-layer GPCANET with varying α (higher α’s for fewer layers, see Appendix. A.7), and select the best α and number of layers using validation set. We report test set performance in Table 1. We do not observe much improvement by GPCANET over other models on smaller datasets CORA, CITESEER, PUBMED, but notable gains on the larger ARXIV and PRODUCTS. As such, GPCANET enables shallow model training via tunable hyperparameter α, achieving comparable or better performance.

heading: 5.3 Q2:UNSUPERVISED VS. SEMI-SUPERVISED GPCA
text: The representations generated by unsupervised GPCA does not use any label information from training data. In this work, we have extended GPCA to (semi-)supervised setting with an additional HP, namely β ∈ [0, 1] that trades-off graph regularization due to the actual input graph edges versus the “ghost” ones added through Y Y T . Overfitting can hurt performance when β is too large or when there is a distribution shift between the training and test sets. For ARXIV and PRODUCTS, we empirically observe that β > 0 always degrades performance, possibly because of the distribution difference between the training and test sets as described in OGB (Hu et al., 2020). Therefore we only study the effect of β on CORA, CITESEER and PUBMED. The pool for β > 0 is {0.1, 0.2}.
CORA CITESEER PUBMED
U-GPCA 81.10 (0.00) 71.80 (0.75) 78.78 (0.36)
S-GPCA (ALL β>0) 81.17 (0.27) 73.20 (0.71) 79.40 (0.69)
S-GPCA β=0.1 81.17 (0.27) 72.07 (0.37) 79.40 (0.69) S-GPCA β=0.2 81.90 (0.00) 73.20 (0.71) 78.73 (0.59)
TABLE 1 BEST 81.48 (0.83) 71.80 (0.75) 78.95 (0.36)

heading: 5.4 Q3: GPCANET-INITIALIZATION FOR GCN
text: Finally, we evaluate the effectiveness of GPCANET-initialization for GCN in terms of performance and robustness under different model sizes, i.e. number of layers L or number of training parameters. For comparison, Xavier initialization (Glorot & Bengio, 2010) is also used to initialize GCN.
DATASET L=2 L=3 L=5 L=10 L=15
CORA XAIVER-INIT 80.62 80.62 79.40 76.37 66.07 CORA GPCANET-INIT 81.67 79.50 80.90 79.82 78.00
CITESEER XAIVER-INIT 71.25 70.15 71.10 61.90 57.40 CITESEER GPCANET-INIT 71.27 69.27 70.15 68.67 67.87
PUBMED XAIVER-INIT 78.42 77.90 77.07 77.00 45.80 PUBMED GPCANET-INIT 78.05 77.25 78.07 77.80 78.03
ARXIV XAIVER-INIT 69.61 70.64 70.33 68.32 61.68 ARXIV GPCANET-INIT 69.76 70.72 70.52 69.77 66.28
PRODUCTS XAIVER-INIT 77.90 78.65 78.08 76.27 74.70 PRODUCTS GPCANET-INIT 78.13 78.71 78.22 77.47 75.90
ing robustness, by reducing performance variation across different seeds. To this end, we first choose the best configuration for each initialization method based on validation performance, and train the GCN model with the chosen configuration using 100 random seeds.
In Figure 2 we present the histogram of test set accuracy over 100 runs with different seeds for ARXIV. (For results on other datasets, see Appendix. A.8.) For both 2-layer and 15-layer GCN, GPCANETinitialization not only outperforms Xavier-initialization w.r.t. average performance, but also in terms of robustness, achieving much lower performance
variation and few bad outliers, especially for deeper GCN. As such, it acts as a strong data-driven prior, facilitating the training of numerous parameters across many layers by identifying a promising region of the parameter space from which supervised fine-tuning is initiated.

heading: 6 CONCLUSION
text: In this work we have (1) discovered a mathematical connection between GPCA and graph convolution of GCN and PPNP; (2) extended GPCA to the (semi-)supervised setting; (3) proposed GPCANET, by stacking GPCA and nonlinear activation, which is a generalized GCN model with an additional hyperparameter to control the degree of graph regularization, and (4) introduced the GPCANETinitialization based on the established connection. Accordingly, we designed extensive experiments demonstrating that (i) the unsupervised shallow GPCA achieves comparable or better performance than GCN, APPNP, as well as other modern GNNs which suggests that graph convolution’s power is mainly driven by graph regularization; (ii) semi-supervised GPCA helps improve performance and should be a powerful yet simple baseline in future research; (iii) GPCANET enables the training of shallow models with competitive performance via increasing the degree of graph regularization at each layer, with reduced memory and training time cost; and finally (iv) GPCANET-initialization acts as a strong data-driven prior for GCN training, enabling robust performance. Our methodological contributions ( 3) & 4) above) capitalize on the discovery of our theoretical findings ( 1) & 2) ), shedding new light toward a better understanding and design of GNNs.

heading: A APPENDIX

heading: A.1 PROOF OF THEOREM 3.1
text: Proof. We give the proof in two steps.
Step 1: For a fixed W , Solve optimal Z∗ as a function of W : When fixing W as constant, the problem becomes quadratic and convex. There is a unique solution, given by first-order optimal condition. Let ℓ denote the objective function as given in equation 5. Its gradient can be calculated as
∂ℓ ∂Z = 2(I + αL̃)Z − 2XW . (15)
Setting equation 15 to 0 leads to the solution Z∗ = (I + αL̃)−1XW .
Step 2: Replace Z with Z∗, Solve optimal W ∗: Substituting Z in objective ℓ with Z∗ = (I + αL̃)−1XW , we reduce the optimization to
min W,WTW=I
∥X− (I+αL̃)−1XWWT ∥2F + αTr [ WTXT (I+αL̃)−1L̃(I+αL̃)−1XW ] . (16)
For this part only, let M = (I + αL̃)−1 to simplify notation. We can show that equation 16 is equivalent to
min W,WTW=I
Tr(XXT +MXWWTWWTXTM)
− 2Tr(MXWWTXT ) + αTr(WTXTML̃MXW ) (17) Using the cyclic property of (Tr)ace (and plugging (I + αL̃)−1 for M back), we can write it as (see Supp. A.2 for detailed derivation.)
max W,WTW=I
Tr [ WTXT (I + αL̃)−1XW ] . (18)
Based on the spectral theorem of PSD matrices, the optimal solution W ∗ of problem equation 18 is the combination of eigenvectors, associated with the largest c eigenvalues of the graph-revised covariance matrix XT (I + αL̃)−1X .

heading: A.2 DERIVATION FROM EQ. EQUATION 17 TO EQ. EQUATION 18
text: For this part only, let A = (I + αL̃)−1 to simplify the notation. We can show that equation 16 is equivalent to
min W,WTW=I
Tr(XXT )− 2Tr(AXWWTXT )
+ Tr(AXWWTWWTXTA) + αTr(WTXTAL̃AXW )
≡ max W,WTW=I 2Tr(AXWWTXT )− Tr(AXWWTXTA)
− αTr(WTXTAL̃AXW ) (19) Using the cyclic property of (Tr)ace, we can write
max W,WTW=I
2Tr(WTXTAXW )− Tr(WTXTAAXW )
− αTr(WTXTAL̃AXW ) max
W,WTW=I Tr
[ WTXT (2A−AA−A(αL̃)A)XW ] max
W,WTW=I Tr
[ WTXT ( A+ {I −A(I + αL̃)}A ) XW ] max
W,WTW=I Tr
[ WTXT (I + αL̃)−1XW ] (20)
where the objective simplifies upon replacing A with (I + αL̃)−1.

heading: A.3 DERIVATION OF EQUIVALENCE IN EQ. EQUATION 9
text: max z
[ corr(Y, z) ]T [ corr(Y, z) ] var(z)
≡ max z
var(Y ) [ corr(Y, z) ]T [ corr(Y, z) ] var(z) (21)
≡ max z
[ cov(Y, z) ]T [ cov(Y, z) ] (22)
where cov(Y, z) = √ var(Y )corr(Y, z) √ var(z)
≡ max z
[ Y T z ]T [ Y T z ] (23)
≡ max z zTY Y T z (24)
Note that in equation 21 we added the term var(Y ) without affecting the optimization problem as it is with respect to z.

heading: A.4 DATASET STATISTICS
text: Datasets used in the experiments are presented in Table 4. Cora, CiteSeer, and PubMed can be downloaded in Pytorch Geometric Library Fey & Lenssen (2019). Arxiv and Products can be accessed in https://ogb.stanford.edu/.

heading: A.5 HYPERPARAMETER CONFIGURATIONS
text: We setup hyperparameters pool for each dataset, presented in Table 5. All methods use the same pool. The only exception is GPCA, as GPCA is just a 1-layer shallow model which can be trained with lager learning rate; we use 0.1 learning rate for it on all datasets.
CORA 0.001 [0.0005, 0.005, 0.05] [2, 3, 5, 10, 15] [128, 256] CITESEER 0.001 [0.0005, 0.005, 0.05] [2, 3, 5, 10, 15] [128, 256] PUBMED 0.001 [0.0005, 0.005, 0.05] [2, 3, 5, 10, 15] [128, 256] ARXIV 0.005 0 [2, 3, 5, 10, 15] [128, 256] PRODUCTS 0.001 0 [2, 3, 5, 10, 15] [128, 256]
DATASET DROPOUT α β
CORA [0, 0.5] [1, 5, 10, 20, 50] [0, 0.1, 0.2] CITESEER [0, 0.5] [1, 5, 10, 20, 50] [0, 0.1, 0.2] PUBMED [0, 0.5] [1, 5, 10, 20, 50] [0, 0.1, 0.2] ARXIV [0, 0.2] [1, 5, 10, 20, 50] 0 PRODUCTS [0, 0.1] [1, 5, 10, 20, 50] 0
Models are trained on every configuration across HP pools and picked based on validation performance. We use the Adam optimizer for all models. Learning rate is first manually tuned for each dataset to achieve stable training, and the same learning rate is fixed for all models—we empirically
observed that learning rate is sensitive to datasets but insensitive to models. For GPCA and GPCANET, number of power iterations in Eq. equation 13 is always set to 5. All experiments use the maximum training epoch as 1000 and repeat 5 times. Detailed configuration of HPs can be found in Supp. A.5. We mainly use a single GTX-1080ti GPU for small datasets CORA, CITESEER, and PUBMED. RTX-3090 GPU is used for ARXIV and PRODUCTS.
Mini-batch training. As nodes are not independent, GNN is mostly trained in full-batch under semi-supervised setting. We use full-batch training for all datasets except PRODUCTS, which is too large to fit into GPU memory during training. ClusterGCN Chiang et al. (2019), a subgraph based mini-batch training algorithm, is used to train GCN and GPCANET. For evaluation, we still use full-batch since a single forward pass can be conducted without memory issues. Initialization is also employed in full-batch.
Fair evaluation. Instead of picking the hyperparameter configurations manually, reported (test) performance is based on the best configuration selected using validation performance, where all models leverage the same hyperparameter pools. Further, each configuration from the pool is conducted 5 times to reduce randomness.
A.6 GPCA WITH VARYING α

heading: A.7 CONFIGURATIONS FOR EXPERIMENTS OF 1∼3-LAYER GPCANET
text: To train a shallow GPCANET with tunable α (β=0 is used), we setup different α pool for different number of layers, because the effect of increasing α is the same to increasing number of layers (shown in Figure 1). We report the pool for α for each layer in Table 7. For other parameters we use the same setting mentioned in Table 5.

heading: A.8 GPCANET-INIT’S ROBUSTNESS FOR ADDITIONAL DATASETS
text: Histogram of test set accuracy over 100 runs for GCN initialized by Xavier-initialization and GPCANET-initialization in CORA (Figure 3), CITESEER (Figure 4), and PUBMED (Figure 5). We have ignored PRODUCTS as it takes too long to run 100 times, but the result should be similar.

heading: A.9 TRAINING CURVE COMPARISON FOR GPCANET-INIT

heading: A.10 RUNTIME COMPARISON
text: We have analyzed the runtime complexity of GPCA, GPCANET, and GPCANET-Init in Sec.3.6 and show their runtime is linear in number of nodes. Table 8 presents the practical runtime comparison among all methods, measured in seconds/epoch for all models, and total initialization seconds for GPCANET-Init, which verified the complexity analysis. Besides, GPCA is a extremely fast method with strong performance, and should be used as a strong baseline in future research.

Review summaries:
This paper relates GCN to PCA from the perspective of optimization. The authors propose Graph PCA that is a general form of GCN. They further introduce a regularization term that enforces nodes with same labels close to each other. 

This manuscript looks at the classic graph-regularized PCA (GPCA) and try to build the connection with GPCA and the state-of-the-art GCN, and finally proposes a new deep graph network GPCANet.  The authors make a number of clear contributions as listed in the paper: 1) they build the connection between GPCA and GCN, 2) Based on this connection, they propose novel way of using such GPCA as a graph layer or as an initialization process for training GCN etc. and 3) thus present a new architecture of GPCANet that performs well on their tests.