Year: 2019

Abstract:
Graph convolutional networks (GCNs) are powerful tools for graph-structured data. However, they have been recently shown to be vulnerable to topological attacks. To enhance adversarial robustness, we go beyond spectral graph theory to robust graph theory. By challenging the classical graph Laplacian, we propose a new convolution operator that is provably robust in the spectral domain and is incorporated in the GCN architecture to improve expressivity and interpretability. By extending the original graph to a sequence of graphs, we also propose a robust training paradigm that encourages transferability across graphs that span a range of spatial and spectral characteristics. The proposed approaches are demonstrated in extensive experiments to simultaneously improve performance in both benign and adversarial situations.

heading: 1 INTRODUCTION
text: Graph convolutional networks (GCNs) are powerful extensions of convolutional neural networks (CNN) to graph-structured data. Recently, GCNs and variants have been applied to a wide range of domains, achieving state-of-the-art performances in social networks (Kipf & Welling, 2017), traffic prediction (Rahimi et al., 2018), recommendation systems (Ying et al., 2018), applied chemistry and biology (Kearnes et al., 2016; Fout et al., 2017), and natural language processing (Atwood & Towsley, 2016; Hamilton et al., 2017; Bastings et al., 2017; Marcheggiani & Titov, 2017), just to name a few (Zhou et al., 2018; Wu et al., 2019).
GCNs belong to a family of spectral methods that deal with spectral representations of graphs (Zhou et al., 2018; Wu et al., 2019). A fundamental ingredient of GCNs is the graph convolution operation defined by the graph Laplacian in the Fourier domain:
gθ ? x := ĝθ(L)x, (1) where x ∈ Rn is the graph signal on the set of vertices V and ĝθ is a spectral function applied to the graph Laplacian L := D −A (whereD andA are the degree matrix and the adjacency matrix, respectively). Because this operation is computational intensive for large graphs and non-spatially localized (Bruna et al., 2014), early attempts relied on a parameterization with smooth coefficients (Henaff et al., 2015) or a truncated expansion in terms of of Chebyshev polynomials (Hammond et al., 2011). By further restricting the Chebyshev polynomial order by 2, the approach in (Kipf & Welling, 2017) referred henceforth as the vanilla GCN pushed the state-of-the-art performance of semi-supervised learning. The network has the following layer-wise update rule:
H(l+1) := ψ ( AH(l)W (l) ) , (2)
where H(l) is the l-th layer hidden state (with H(1) := X as nodal features), W (l) is the l-th layer weight matrix, ψ is the usual point-wise activation function, and A is the convolution operator chosen to be the degree weighted Laplacian with some slight modifications (Kipf & Welling, 2017). Subsequent GCN variants have different architectures, but they all share the use of the Laplacian matrix as the convolution operator (Zhou et al., 2018; Wu et al., 2019).

heading: 1.1 WHY NOT GRAPH LAPLACIAN?
text: Undoubtedly, the Laplacian operator (and its variants, e.g., normalized/powered Laplacian) plays a central role in spectral theory, and is a natural choice for a variety of spectral algorithms such as
principal component analysis, clustering and linear embeddings (Chung & Graham, 1997; Belkin & Niyogi, 2002). So what can be problematic?
From a spatial perspective, GCNs with d layers cannot acquire nodal information beyond its ddistance neighbors; hence, it severely limits its scope of data fusion. Recent works (Lee et al., 2018; Abu-El-Haija et al., 2018; 2019; Wu et al., 2019) alleviated this issue by directly powering the graph Laplacian.
From a spectral perspective, one could demand better spectral properties, given that GCN is fundamentally a particular (yet effective) approximation of the spectral convolution (1). A key desirable property for generic spectral methods is known as “spectral separation,” namely the spectrum should comprise a few dominant eigenvalues whose associated eigenvectors reveal the sought structure in the graph. A well-known prototype is the Ramanujan property, for which the second leading eigenvalue of a r-regular graph is no larger than 2 √ r − 1, which is also enjoyed asymptotically by random r-regular graphs (Friedman, 2004) and Erdős-Rényi graphs that are not too sparse (Feige & Ofek, 2005). In a more realistic scenario, consider the stochastic block model (SBM), which attempts to capture the essence of many networks, including social graphs, citation graphs, and even brain networks (Holland et al., 1983).
Definition 1 (Simplified stochastic block model). The graph G with n nodes is drawn under SBM(n, k, aintra, ainter) if the nodes are evenly and randomly partitioned into k communities, and nodes i and j are connected with probability aintra/n ∈ [0, 1] if they belong to the same community, and ainter/n ∈ [0, 1] if they are from different communities.
It turns out that for community detection, the top k leading eigenvectors of the adjacency matrix A play an important role. In particular, for the case of 2 communities, spectral bisection algorithms simply take the second eigenvector to reveal the community structure. This can be also seen from the expected adjacency matrix E[A] under SBM(n, 2, aintra, ainter), which is a rank-2 matrix with the top eigenvalue 12 (aintra + ainter) and eigenvector 1, and the second eigenvalue 1 2 (aintra − ainter) and eigenvector σ such that σi = 1 if i is in community 1 and σi = −1 otherwise. More generally, the second eigenvalue is of particular theoretical interests because it controls at the first order how fast heat diffuses through graph, as depicted by the discrete Cheeger inequality (Lee et al., 2014).
While one would expect taking the second eigenvector of the adjacency matrix suffices, it often fails in practice (even when it is theoretically possible to recover the clusters given the signal-to-noise ratio). This is especially true for sparse networks, whose average nodal degrees is a constant that does not grow with the network size. This is because the spectrum of the Laplacian or adjacency matrix is blurred by “outlier” eigenvalues in the sparse regime, which is often caused by high degree nodes (Kaufmann et al., 2016). Unsurprisingly, powering the Laplacian would be of no avail, because it does not change the eigenvectors or the ordering of eigenvalues. In fact, those outliers can become more salient after powering, thereby weakening the useful spectral signal even further. Besides, pruning the largest degree nodes in the adjacency matrix or normalizing the Laplacian cannot solve the issue. To date, the best results for pruning does not apply down to the theoretical recovery threshold (Coja-Oghlan, 2010; Mossel et al., 2012; Le et al., 2015); either outliers would persist or one could prune too much that the graph is destroyed. As for normalized Laplacian, it may overcorrect the large degree nodes, such that the leading eigenvectors would catch the “tails” of the graph, i.e., components weakly connected to the main graph. See Figure A.3 in the appendix for an experimental illustration.

heading: In summary, graph Laplacian may not be the ideal choice due to its limited spatial scope of information fusion, and its undesirable artefacts in the spectral domain.

heading: 1.2 IF NOT LAPLACIAN, THEN WHAT?
text: In searching for alternatives, potential choices are many, so it is necessary to clarify the goals. In view of the aforementioned pitfalls of the graph Laplacian, one would naturally ask the question:
Can we find an operator that has wider spatial scope, more robust spectral properties, and is meanwhile interpretable and can increase the expressive power of GCNs?
From a perspective of graph data analytics, this question gauges how information is propagated and fused on a graph, and how we should interpret “adjacency” in a much broader sense. An image can
be viewed as a regular grid, yet the operation of a CNN filter goes beyond the nearest pixel to a local neighborhood to extract useful features. How to draw an analogy to graphs?
From a perspective of robust learning, this question sheds light on the basic observation that realworld graphs are often noisy and even adversarial. The nice spectral properties of a graph topology can be lost with the presence or absence of edges. What are some principled ways to robustify the convolution operator and graph embeddings?
In this paper, we propose a graph learning paradigm that aims at achieving this goal, as illustrated in Figure 1. The key idea is to generate a sequence of graphs from the given graph that capture a wide range of spectral and spatial behaviors. We propose a new operator based on this derived sequence. Definition 2 (Variable power operator). Consider an unweighted and undirected graph G. Let A[k] denote the k-distance adjacency matrix, i.e., [ A[k] ] ij
= 1 if and only if the shortest distance (in the original graph) between nodes i and j is k. The variable power operator of order r is defined as:
A (r) θ = r∑ k=0 θkA [k], (3)
where θ := (θ0, . . . , θr) is a set of parameters.
Clearly, A(r)θ is a natural extension of the classical adjacency matrix (i.e., r = 1 and θ0 = θ1 = 1). With power order r > 1, one can increase the spatial scope of information fusion on the graph when applying the convolution operation. The parameters θk also has a natural explanation—the magnitude and the sign of θk can be viewed as “global influence strength” and “global influence propensity” at distance k, respectively, which also determines the participation factor of each graph in the sequence in the aggregated operator.
Furthermore, we provide some theoretical justification of the proposed operator by establishing the following asymptotic property of spectral separation under the important SBM setting, which is, nevertheless, not enjoyed by the classical Laplacian operator or its normalized or powered versions. (All proofs are given in the appendix.) Theorem 3 (Asymptotic spectral separation of variable power operator). Consider a graph G drawn from SBM(n, 2, aintra, ainter). Assume that the signal-to-noise ratio ξ22/ξ1 > 1, where ξ1 = 1 2 (aintra + ainter) and ξ2 = 12 (aintra − ainter) (c.f., (Decelle et al., 2011)). Suppose r is on the order of c log(n) for a constant c, such that c log(ξ1) < 1/4. Given nonvanishing θk for k > r/2, the variable power operator A(r)θ has the following spectral properties: (i) the leading eigenvalue is on the order of Θ (‖θ‖1ξr1), the second leading eigenvalue is on the order of Θ (‖θ‖1ξr2), and the rest are bounded by ‖θ‖1n ξr/21 O(log(n)) for any fixed > 0; and (ii) the two leading eigenvectors are sufficient to recover the two communities asymptotically (i.e., as n goes to infinity).
Intuitively, the above theoretical result suggests that the variable power operator is able to “magnify” benign signals from the latent graph structure while “suppressing” noises due to random artefacts.
This is expected to improve spectral methods in general, especially when the benign signals tend to be overwhelmed by noises. For the rest of the paper, we will apply this insight to propose a robust graph learning paradigm in Section 2, as well as a new GCN architecture in Section 3. We also provide empirical evidence of the gain from this theory in Section 4 and conclude in Section 5.

heading: 1.3 RELATED WORK
text: Beyond nearest neighbors. Several works have been proposed to address the issue of limited spatial scope by powering the adjacency matrix (Lee et al., 2018; Wu et al., 2019; Li et al., 2019). However, simply powering the adjacency does not extract spectral gap and may even make the eigenspectrum more sensitive to perturbations. Abu-El-Haija et al. (2018; 2019) also introduced weight matrices for neighbors at different distances. But this could substantially increase the risk of overfitting in the low-data regime and make the network vulnerable to adversarial attacks.
Robust spectral theory. The robustness of spectral methods has been extensively studied for graph partitioning/clustering (Li et al., 2007; Balakrishnan et al., 2011; Chaudhuri et al.; Amini et al., 2013; Joseph et al., 2016; Diakonikolas et al., 2019). Most recently, operators based on self-avoiding or nonbacktracking walks have become popular for SBM (Massoulié, 2014; Mossel et al., 2013; Bordenave et al., 2015), which provably achieve the detection threshold conjectured by Decelle et al. (2011). Our work is partly motivated by the graph powering approach by Abbe et al. (2018), which leveraged the result of Massoulié (2014); Bordenave et al. (2015) to prove the spectral gap. The main difference with this line of work is that these operators are studied only for spectral clustering without incorporating nodal features. Our proposed variable power operator can be viewed as a kernel version of the graph powering operator (Abbe et al., 2018), thus substantially increasing the capability of learning complex nodal interactions while maintaining the spectral property.
Robust graph neural network. While there is a surge of adversarial attacks on graph neural networks (GNNs) (Dai et al., 2018; Zügner & Günnemann, 2019; Bojchevski & Günnemann, 2019), very few methods have been proposed for defense (Sun et al., 2018). Existing works employed known techniques from computer vision (Szegedy et al., 2014; Goodfellow et al., 2015; Szegedy et al., 2016), such as adversarial training with “soft labels” (Chen et al., 2019) or outlier detection in the hidden layers (Zhu et al., 2019), but they do not exploit the unique characteristics of graphstructured data. Importantly, our approach simultaneously improves performance in both the benign and adversarial tests, as shown in Figure 2 (details are presented in Section 4).

heading: 2 GRAPH AUGMENTATION: ROBUST TRAINING PARADIGM
text: Exploration of the spectrum of spectral and spatial behaviors. Given a graph G, consider a family of its powered graphs, { G(2), . . . ,G(r) } , where G(k) is obtained by connecting nodes with distance less than or equal to k. This “graph augmentation” procedure is similar to “data augmentation”, because instead of limiting the learning on a single graph that is given, we artificially create a series of graphs that are closely related to each other in the spatial and spectral domains.
As we increase the power order, the graph also becomes more homogenized (Figure A.1). In particular, it can help near-isolated nodes (i.e., low-degree vertices), since they become connected beyond their nearest neighbors. By comparison, simply raising the adjacency matrix to its r-th power will
make them appear even more isolated, because it inadvertently promotes nodes with high degrees or nearby loops much more substantially as a result of feedback. Furthermore, the powered graphs can extract spectral gaps in the original graph despite local irregularities, thus boosting the signalto-noise ratio in the spectral domain.
Transfer of knowledge from the powered graph sequence. Consider a generic learning task on a graph G with data D. The loss function is denoted by `(W;G,D) for a particular GCN architecture parametrized byW . For instance, in semi-supervised learning,D consists of features and labels on a small set of nodes (see Table A.1), and ` is the cross-entropy loss over all labeled examples. Instead of minimizing over `(W;G,D) alone, we use all the powered graphs:
`(W;G,D) + r∑
k=2
αk`(W;G(k),D), (r-GCN)
where αk ≥ 0 gauges how much information one desires to transfer from powered graph G(k) to the learning process. By minimizing the (r-GCN) objective, one seeks to optimize the network parameterW on multiple graphs simultaneously, which is beneficial in two ways: (i) in the low-data regime, like semi-supervised learning, it helps to reduce the variance to improve generalization and transferability; (ii) in the adversarial setting, it robustifies the network since it is more likely that the perturbed network is contained in the wider spectrum during training.

heading: 3 FROM FIXED TO VARIABLE POWER NETWORK
text: By using the variable power operator illustrated in Figure 1, we substantially increase the search space of graph operators. The proposed operator also leads to a new family of graph algorithms with broader scope of spatial fusion and enhanced spectral robustness. As the power grows, the network eventually becomes dense. To manage this increased complexity and make the network more robust against adversarial attacks in the feature domain, we propose a pruning mechanism.
Graph sparsification. Given a graph G := (V, E [1]), consider its powered version G(r) := (V, E(r)) and a sequence of intermediate graphs G[2], . . . ,G[r], where G[k] := (V, E [k]) is constructed by connecting two vertices if and only if the shortest distance is k in G. Clearly, {E [k]}rk=1 forms a partition of E(r). For each node i ∈ V , denote its r-neighborhood by Nr(i) := {j ∈ V | dG(i, j) ≤ r}, which is identical to the set of nodes adjacent to i in G[r]. Next, for each edge within this neighborhood, we associate a value using some suitable distance metric φ to measure “aloofness.” For instance, it can be the usual Euclidean distance or correlation distance in the feature space. Based on this formulation, we prune an edge e := (i, j) in E(r) if the value is larger than a threshold τ(i, j), and denote the edge set after pruning Ē(r). Then, we can construct a new sequence of sparsified graphs, G[k] with edge sets Ē [k] = E [k] ∩ Ē(r) and adjacency matrix Ā[k]. Hence, the variable power operator is given by Ā(r)θ = ∑r k=0 θkĀ
[k]. Due to the influence of high-degree nodes in the spectral domain, one can adaptively choose the thresholds τ(i, j) to mitigate their effects. Specifically, we choose τ to be a small number if either i or j are nodes with high degrees, thereby making the sparsification more influential in densely connected neighborhoods than weakly connected parts.
Layer-wise update rule. To demonstrate the effectiveness of the proposed operator, we adopt the vanilla GCN strategy (2). Importantly, we replace the graph convolutional operator A with the variable power operator to obtain the variable power network (VPN):
A = D− 1 2 (I + Ā
(r) θ )D
− 12 , (VPN)
where Dii = 1 + |{j ∈ V | dG(i, j) = 1}|. The introduction of I is reminiscent of the “renormalization trick” (Kipf & Welling, 2017), but it can be also viewed as a regularization strategy in this context, which is well-known to improve the spectral robustness (Amini et al., 2013; Joseph et al., 2016). This construction immediately increases the scope of data fusion by a factor of r. Proposition 4. By choosingA with (VPN) in the layer-wise update rule (2), the output at each node from a L-layer GCN depends on neighbors within L ∗ r hops.
Since we proved that the variable power operator has nice spectral separation in Theorem 3, VPN is expected to promote useful spectral signals from the graph topology (similar to the preservation of
useful information in images (Jacobsen et al., 2018), our method preserves useful information in the graph topology). This claim is substantiated with the following proposition. Proposition 5. Given a graph with two hidden communities. Consider a 2-layer GCN architecture with layer-wise update rule (2). Suppose that A has a spectral gap. Further, assume that the leading two eigenvectors are asymptotically aligned with 1 and ν, i.e., the community membership vector, and that both are in the range of feature matrix X . Then, there exists a configuration of W (1) and W (2) such that the GCN outputs can recover the community with high probability.

heading: 4 EXPERIMENTS
text: The proposed methods are evaluated against several recent models, including vanilla GCN (Kipf & Welling, 2017) and its variant PowerLaplacian where we simply replace the adjacent matrix with its powered version, three baselines using powered Laplacian IGCN (Li et al., 2019), SGC (Wu et al., 2019) and LNet (Liao et al., 2019), the recent method MixHop (Abu-El-Haija et al., 2019) which attempts to increase spatial scope fusion, as well as a state-of-the-art baseline and RGCN (Zhu et al., 2019), which is also aimed at improving the robustness of Vanilla GCN. All baseline methods on based on their public codes.

heading: 4.1 REVISITING THE STOCHASTIC BLOCK MODEL
text: SBM dataset. We generated a set of networks under SBM with 4000 nodes and parameters such that the SNRs range from 0.58 to 0.91. To disentangle the effects from nodal features with that from the spectral signal, we set the nodal features to be one-hot vectors. The label rates are 0.1%, 0.5% and 1%, the validation rate is 1%, and the rest of the nodes are testing points.
Performance. Since the nodal features do not contain any useful information, learning without topology such as multi-layer perceptron (MLP) is only as good as random guessing. The incorporation of graph topology improves classification performance—the higher the SNR (i.e., ξ22/ξ1, see Theorem 3), the higher the accuracy. Overall, as shown in Figure 3, the performance of the proposed method (VPN) is superior than other baselines, which either use Laplacian (GCN, Chebyshev, RGCN, LNet) or its powered variants (IGCN, PowerLaplacian).
Spectral separation and Fourier modes. From the eigenspectrum of the convolution operators (Figure A.3), we see that the spectral separation property is uniquely possessed by VPN, whose first two leading eigenvectors carry useful information about the underlying communities: without the help of nodal features, the accuracy is 87% even with label rate of 0.5%. Let Φ denote the Fourier
modes of the the adjacency matrix or VPN, and X be the nodal features (i.e., identity matrix). We analyze the information from spectral signals (e.g., the k-th and k+1-th eigenvectors) by estimating the accuracy of an MLP with filtered nodal features, namely Φ:,k:(k+1)Φ>:,k:(k+1)X , as shown in Figure 4. The accuracy reflects the information content in the frequency components. We see that the two leading eigenvectors of VPN are sufficient to perform classification, whereas those of the adjacency matrix cannot make accurate inferences.

heading: 4.2 SEMI-SUPERVISED NODE CLASSIFICATION
text: Experimental setup. We followed the setup of (Yang et al., 2016; Kipf & Welling, 2017) for citation networks Citeseer, Cora and Pubmed (please refer to the Appendix for more details).
Graph powering order can influence spatial and spectral behaviors. Our theory suggests powering to the order of log(n); in practice, orders of 2 to 4 suffice (Figure 5). Here, we chose the power order to be 4 for r-GCN on Citeseer and Cora, and 3 for Pubmed, and reduced the order by 1 for VPN.
Performance. By replacing Laplacian with VPN, we see an immediate improvement in performance (Table 1). We also see that a succinct parametrization of the global in-
fluence relation in VPN is able to increase the expressivity of the network. For instance, the learned θ at distances 2 and 3 for Citeseer are 3.15e-3 and 3.11e-3 with p-value less than 1e-5. This implies that the network tends to put more weights in closer neighbors.

heading: 4.3 DEFENSE AGAINST EVASION ATTACKS
text: To evaluate the robustness of the learned network, we considered the setting of evasion attacks, where the model is trained on benign data but tested on adversarial data.
Adversarial methods. Five strong global attack methods are considered, including DICE (Zügner & Günnemann, 2019), Aabr and ADW3 (Bojchevski & Günnemann, 2019), Meta-Train and MetaSelf (Zügner & Günnemann, 2019). We further modulated the severity of attack methods by varying the attack rate, which corresponds to the percentage of edges to be attacked.
Robustness evaluation. In general, both r-GCN and VPN are able to improve over baselines for the defense against evasion attacks, e.g., Figure 6 for theADW3 attack (detailed results for other attacks are listed in the Appendix). It can be also observed that the proposed methods are more robust in Citeseer and Cora than Pubmed. In addition to the low label rates, we conjectured that topological
Attack rate (%)
R ob
us tn es s
en ha
nc em
en t( % )
Attack rate (%)
Figure 6: Robustness merit
(
accuracypost-attackproposed-method − accuracy post-attack vanilla GCN
) reported in percentage un-
der ADW3 attack. The error bar indicates standard deviation over 20 independent simulations.
attacks are more difficult to defend for networks with prevalent high-degree nodes, because the attacker can bring in more irrelevant vertices by simply adding a link to the high-degree nodes.
Informative and robust low-frequency spectral signal. It has been observed by Wu et al. (2019); Maehara (2019) that GCNs share similar characteristics of low-pass filters, in the sense that nodal features filtered by low-frequency Fourier modes lead to accurate classifiers (e.g., MLP). However, one key question left unanswered is how to obtain the Fourier modes. In their experiments, they derive it from the graph Laplacian. By using VPN to construct the Fourier modes, we show that the information content in the low-frequency domain can be improved.
More specifically, we first perform eigendecomposition of the graph convolutional operator (i.e., graph Laplacian or VPN) to obtain the Fourier modes Φ. We then reconstruct the nodal features X using only the k-th and the k + 1-th eigenvectors, i.e., Φ:,k:(k+1)Φ>:,k:(k+1)X . We then use the reconstructed features in MLP to perform the classification task in a supervised learning setting. As Figure 7 shows, features filtered by the leading eigenvectors of VPN lead to higher classification accuracy compared to the classical adjacency matrix.
For the adversarial testing, we construct a new basis Φ̃ based on the attacked graph, and then use Φ̃:,k:(k+1)Φ̃ > :,k:(k+1)X as test points for the MLP trained in the clean data setting. As can be seen in Figure 8, models trained based on VPN filtered features also have better adversarial robustness in evasion attacks. Since the eigenvalues of the corresponding operator exhibit low-pass filtering characteristics (Figure A.2), the enhanced benign and adversarial accuracy of VPN is attributed to the increased signal-to-noise ratio in the low-frequency domain. This is in alignment with the theoretical proof of spectral gap developed in this study.

heading: 4.4 DISCUSSIONS
text: In this paper, we proposed a new graph operator to replace the graph Laplacian in the classic GCN. It is worth mentioning that our approach is very different from the k-th order polynomials of the Laplacian developed in (Defferrard et al., 2016) or the high-order graph Laplacian employed in (Lee et al., 2018; Li et al., 2019; Abu-El-Haija et al., 2019; Wu et al., 2019). If the graph Laplacian L is powered to the k-th order, the resulting matrix has the same eigenvectors asL, with only the eigenvalues
Φ̃:,1:kΦ̃ > :,1:kX as reconstructed nodal features for testing, where k is 10 for Citeseer and 5 for Cora and Pubmed. Error bar indicates standard deviation over 10 runs.
powered to the k-th order. Thus, the corresponding k-th order polynomial has the same eigenvectors, with the polynomial function applied only to the eigenvalues. However, the proposed variable power operator has radically different eigenvectors as the graph Laplacian or its powered variant. This is an important difference because we observed empirically that the leading eigenvectors of a graph Laplacian are extremely sensitive to outliers under the SBM, and they often correspond to either tails or high-degree nodes (please see Figure A.3 for an illustration). However, our result in Theorem 3 suggests that the leading eigenvectors of the proposed operator asymptotically recovers the underlying community under SBM, and enjoys the “spectral gap” property.
We also remark that the meaning of robustness in the machine learning literature (“adversarial robustness”) and in the spectral theory literature (“spectral robustness”) appear to be different. However, in this work, we argue that spectral robustness is closely related to adversarial robustness. This can be explained from a matrix perturbation point of view. Let λk and φk be the eigenvalue and the corresponding eigenvector of the adjacency matrixA. Then, the sensitivity of the eigenvector to the perturbation is given by (Trefethen & Bau III, 1997):
∂φk ∂A(ij) = n∑ `=1, 6̀=k φ`(i)φk(j)(2− δ(ij)) λk − λ` φ`,
where A(ij) is the (i, j)-th entry of A, φk(j) is the j-th entry of φk, δ(ij) is the perturbation of the (i, j)-th entry of A, and n is the dimension of A. It can be seen that the perturbation of the convolutional operator is controlled by the inverse of the spectral gap. When the spectral gaps between the leading eigenvalues (i.e., λ1 and λ2) and the rest of the eigenvalues are large, the perturbations of the leading eigenvectors due to the adversarial attack are controlled. Since the eigenvalues of the graph convolutional operator decays very quickly (see Figure A.2 as an illustration), the learned representation is heavily influenced by the leading eigenvectors. Thus, by controlling the perturbations of the leading eigenvectors, one can expect to control the perturbations of the outputs.
Our theoretical result in Theorem 3 was based on the standard SBM, which is a classic model proposed in mathematical sociology and has been adopted to model and analyze real-world social and biological networks (Funke & Becker, 2019). The main advantage of SBM is that it encodes the structural information that nodes belong to the same community tend to be more connected with each other. Nevertheless, we recognize one limitation of our analysis is that SBM might not be a suitable model for some existing or future applications, in which case, models with more flexibility in the degree distributions, such as the degree-corrected SBM proposed by Karrer & Newman (2011), labeled SBM by Heimlicher et al. (2012), or hiearchical SBM proposed by Peixoto (2017) could be more applicable. Also our analysis is limited to two communities, and the extension to more than two communities is a challenging open problem. Nevertheless, through extensive experiments, we observed that our theory has strong implications for real-world graphs. This is evident from the informative and robust low-frequency spectral signal results in Figures 7 and 8, which show that the proposed graph convolution operator has a clear advantage in the low-frequency regime over the classic graph Laplacian in terms of classification accuracy.

heading: 5 CONCLUSION
text: This study goes beyond classical spectral graph theory to defend GCNs against adversarial attacks. We challenge the central building block of existing methods, namely the graph Laplacian, which
is not robust to noisy links. For adversarial robustness, spectral separation is a desirable property. We propose a new operator that enjoys this property and can be incorporated in GCNs to improve expressivity and interpretability. Furthermore, by generating a sequence of powered graphs based on the original graph, we can explore a spectrum of spectral and spatial behaviors and encourage transferability across graphs. The proposed methods are shown to improve both benign and adversarial accuracy over various baselines evaluated against a comprehensive set of attack strategies.

heading: A APPENDIX
text: A.1 PROOF OF PROPOSITION 4
To validate the claim, we prove a more general result. Consider a graph G and its powered graph G(r), whose adjacency matrices areA andA(r), respectively. Since the support ofA(r)θ is identical to that of A(r), we will focus on A(r) for notational simplicity. For any signal f : V 7→ R|V| on a graph, denote the support of the signal fv on node v ∈ V by Sfv , i.e., fv only depends on the signals on nodes in Sfv . Then, it is easy to see that the signal g , Af has support Sgv ⊆ ∪u∼rvSfu , where u ∼r v denotes adjacency relation on G(r) defined by A(r). Clearly, the set defined by {u ∈ V | u ∼r v} is identical to {u ∈ V | dG(u, v) ≤ r}, where dG(u, v) denotes the shortest distance from u to v on G. Also, we see that element-wise operation does not expand the support, i.e., for any element-wise activation σ, the signal gσ , σ(f) has support Sgσv ⊆ Sfv . In particular, we can prove the claim by iteratively applying the function composition layer by layer, since each layer is of the form σ(A
(r) θ H (l)W (l)), where H(l) is the hidden state at layer l and W (l) is the corresponding weight matrix. Hence after L iterations of layer-wise updates, the support of the output on v is given by {u ∈ V | dG(u, v) ≤ r ∗ L}.
A.2 PROOF OF PROPOSITION 5
Let λ1 and λ2 be the two leading eigenvalues of Ā with corresponding eigenvectors φ1 and φ2, respectively. Without loss of generality, assume that both eigenvalues are nonnegative. Under the
assumption thatφ1 andφ2 lie in the range ofX , there exists aW (1) such thatXW (1) = [φ1 φ2]. So the output of the first layer with ReLU activation is given by:
H(1) = σ ( ĀXW (1) ) = σ ( Ā [φ1 φ2] ) (4)
= [λ1σ(φ1) λ2σ(φ2)] (5) = [ λ1(1 + σ(h1)) λ2( 1 2ν + 1 21 + σ(h2)) ] where ‖h1‖ = o(1) and ‖h2‖ = o(1) (6)
where (5) is due to the definition of eigenvectors and that σ is an element-wise operation, and (6) is due to the asymptotic alignment condition and the fact that σ(ν) = 12ν + 1 21. Therefore, by
choosing the last layer weightW (2) =
[ − 1λ1λ2 1 λ1λ2
2 λ22
− 2 λ22
] , the output of the last layer is given by:
Y = ĀH(1)W (2) = Ā [ λ1(1 + σ(h1)) λ2( 1 2ν + 1 21 + σ(h2)) ] W (2)
= [ λ211 + h ′ 1 λ22 2 ν + λ1λ2 2 1 + h ′ 2 ] W (2)
= [ν + h′′1 −ν + h′′2 ] where ‖h′k‖2 = o(1), ‖h′′k‖2 = o(1) for k = 1, 2. Therefore, after the softmax operation, we can see that the 2-layer output recovers the membership of the elements exactly.
One general remark is that while the above argument is not limited to the leading two eigenvalues/eigenvectors, from the form of W (2) whose elements depend inversely on the corresponding eigenvalues, we can expect more numerical stability (and hence easier learning) when the corresponding eigenvalues are large.
A.3 PROOF OF THEOREM 3
Consider a stochastic block model with two communities, where the intra-community connection probability is aintra/n and the inter-community connection probability is ainter/n, and assume that aintra > ainter. Let ξ1 = aintra+ainter2 and ξ2 = aintra−ainter 2 , which correspond to the first and second eigenvalues of the expected adjacency matrix of the above model. Let ν denote the community label vector, with νi = ±1 depending on the membership of node i. For weak recovery, we assume the detectability condition (Decelle et al., 2011):
ξ22 > ξ1. (7)
With the notations from the main paper, let A(r)θ = ∑r k=0 θkA [k] denote the θ-weighted variable
power operator, where [ A[k] ] ij = { 1 if dG(i, j) = k 0 o.w. is the k-distance adjacency matrix. Let A{k} denote the nonbacktracking path counting matrix, where the (i, j)-th component indicate the number of self-avoiding paths of graph edges of length k connecting i to j (Massoulié, 2014). Our goal is to prove that the top eigenvectors of A(r)θ are asymptotically aligned with those of A
{k} for k greater than log(n) up to a constant multiplicative factor. To do so, we first recall the result of (Massoulié, 2014), which examines the behaviors of top eigenvectors ofA{k}. Lemma 6 ((Massoulié, 2014)). Assume that (7) holds. Then, w.h.p., for all k ∈ {r/2, ..., r}:
(a) The operator norm ‖A{k}‖2 is up to logarithmic factors Θ(ξk1 ). The second eigenvalue of A{k} is up to logarithmic factors Ω(ξk2 ).
(b) The leading eigenvector is asymptotically aligned withA{k}1: A{k}A{k}1
‖A{k}A{k}1‖2 =
A{k}1
‖A{k}1‖2 + hk
where ‖hk‖ = o(1). The second eigenvector is asymptotically aligned withA{k}ν: A{k}A{k}ν
‖A{k}A{k}ν‖2 =
A{k}ν
‖A{k}ν‖2 + h′k
where ‖h′k‖ = o(1).
Now, we first define Γ{r}θ = ∑r k=r/2 θkA {k}, so
A (r) θ = Γ {r} θ + ∆ {r} θ = Γ {r} θ +A (r/2−1) θ + r∑ k=r/2 θk(A [k] −A{k}).
Also, by (Massoulié, 2014, Theorem 2.3), the local neighborhoods of vertices do not grow fast, and w.h.p., the graph of A(r/2−1) has maximum degree O(ξr/21 (log n)
2). Let φ be the leading eigenvector of A(r/2−1)θ , and let v be the node where φ is maximum (i.e., φ(v) ≥ φ(u)). Without loss of generality, assume that φ(v) > 0. We have a good control over the first eigenvalue λ1 of A
(r/2−1) θ as follows:
λ1 = (A
(r/2−1) θ φ)(v)
φ(v) =
∑r/2−1 k=0 ∑ {u:dG(u,v)=k} θkφ(u)
φ(v)
= r/2−1∑ k=0 θk ∑
{u:dG(u,v)=k}
φ(u)
φ(v)
≤ r/2−1∑ k=0 |θk| ∑
{u:dG(u,v)=k}
1
≤ ‖θ‖∞ ∑
{u:dG(u,v)≤r/2−1}
1
= O(‖θ‖1ξr/21 (log n)2)
Due to (Abbe et al., 2018, Theorem 2), which proves an identical result of Lemma 6 but for A[k], we know that w.h.p., for r = log(n), we can bound the last term by
‖ r∑
k=r/2
θk(A [k] −A{k})‖2 ≤ r∑ k=0 |θk|‖A[k] −A{k}‖2 = O ( ‖θ‖1ξr/21 (log n)4 )
Therefore, by triangle inequality, we can conclude that
‖∆{r}θ ‖2 = ‖A (r) θ − Γ {r} θ ‖2 = O ( ‖θ‖1ξr/21 (log n)4 ) . (8)
Now, our goal is to show a similar result as Lemma 6 for Γ{r}θ . In particular,
‖Γ{r}θ A{r}1‖2 ‖A{r}1‖2 = ∥∥∥∥∥∥ r∑
k=r/2
θkA {k} A
{r}1
‖A{r}1‖2 ∥∥∥∥∥∥ 2
(9)
= ∥∥∥∥∥∥ r∑
k=r/2
θkA {k} A
{k}1
‖A{k}1‖2 ∥∥∥∥∥∥ 2 +O  r∑ k=r/2 |θk|‖A{k}‖2 ∥∥∥∥ A{r}1‖A{r}1‖2 − A {k}1 ‖A{k}1‖2 ∥∥∥∥  (10)
= ∥∥∥∥∥∥ r∑
k=r/2
θkA {k} A
{k}1
‖A{k}1‖2 ∥∥∥∥∥∥ 2 +O  r∑ k=r/2 |θk|ξk1n−δ  (11)
= ∥∥∥∥∥∥ r∑
k=r/2
θkA {k} A
{k}1
‖A{k}1‖2 ∥∥∥∥∥∥ 2 +O ( ‖θ‖1ξr1n−δ ) (12)
= ∥∥∥∥∥∥ r∑
k=r/2
θk
( A{k}1
‖A{k}1‖2 + hk ) ‖A{k}A{k}1‖2 ‖A{k}1‖2 ∥∥∥∥∥∥ 2 +O ( ‖θ‖1ξr1n−δ ) for ‖hk‖2 = o(1)
(13)
= ∥∥∥∥∥∥ r∑
k=r/2
θk
( A{r}1
‖A{r}1‖2 + hk ) ‖A{k}A{k}1‖2 ‖A{k}1‖2 ∥∥∥∥∥∥ 2 +
O ∥∥∥∥∥∥ r∑
k=r/2
θk ‖A{k}A{k}1‖2 ‖A{k}1‖2 ∥∥∥∥∥∥ 2 ∥∥∥∥ A{r}1‖A{r}1‖2 − A {k}1 ‖A{k}1‖2 ∥∥∥∥ +O (‖θ‖1ξr1n−δ)
(14)
= ∥∥∥∥∥∥ r∑
k=r/2
θk
( A{r}1
‖A{r}1‖2 + hk ) ‖A{k}A{k}1‖2 ‖A{k}1‖2 ∥∥∥∥∥∥ 2 +O ( ‖θ‖1ξr1n−δ ) (15)
=  r∑ k=r/2 |θk| ‖A{k}A{k}1‖2 ‖A{k}1‖2  (1 + o(1)) +O (‖θ‖1ξr1n−δ) (16) = O ( ‖θ‖1ξr1n−δ ) , (17)
where (10) and (14) are due to triangle inequalities, (11) and (15) are due to results in (Massoulié, 2014), and (13) is due to Lemma 6. Therefore, we have ‖Γ{r}θ A{r}1‖2 = Θ ( ‖θ‖1ξr1‖A{r}1‖2 ) .
Similarly, we can prove that ‖Γ{r}θ A{r}ν‖2 = Θ ( ‖θ‖1ξr2‖A{r}ν‖2 ) .
Furthermore, we can show that Γ{r}θ satisfies the weak Ramanujan property (the proof is similar to (Massoulié, 2014) and is deferred to Section A.4).
Lemma 7. For any fixed > 0, Γ{r}θ satisfies the following weak Ramanujan property w.h.p.,
sup ‖u‖2=1,u>A{r}1=u>A{r}ν=0
‖Γ{r}θ u‖2 = ‖θ‖1n ξ r/2 1 O(log(n))
By Lemma 7, the leading two eigenvectors of Γ{r}θ will be asymptotically in the span ofA {r}1 and A{r}ν according to the variational definition of eigenvalues. By our previous analysis, this means that the top eigenvalue of Γ{r}θ will be Θ (‖θ‖1ξr1) with eigenvector asymptotically aligned with A{r}1. Since by (Massoulié, 2014, Lemma 4.4), A{r}1 and A{r}ν are asymptotically orthogonal, the second eigenvalue of Γ{r}θ will be Θ (‖θ‖1ξr2) with eigenvector asymptotically aligned with A{r}ν.
Since the perturbation term ‖∆{r}θ ‖2 = o ( λ2(Γ {r} θ ) ) by (8), using Weyl’s inequality (Weyl, 1912),
we can conclude that the leading eigenvalue λ1 ( A (r) θ ) = Θ (‖θ‖1ξr1), the second eigenvalue
λ2
( A
(r) θ ) = Θ (‖θ‖1ξr2), and the rest of the eigenvalues are bounded by ‖θ‖1n ξ r/2 1 O(log(n)).
Moreover, since ‖∆{r}θ ‖2 = o ( max { λ1(Γ {r} θ )− λ2(Γ {r} θ ), λ2(Γ {r} θ )− λ3(Γ {r} θ ) }) , by the
Davis-Kahan Theorem (Davis & Kahan, 1970), the leading two eigenvectors of A(r)θ are asymptotically aligned with A{r}1 and A{r}ν, which is shown to be enough for the rounding procedure of (Massoulié, 2014) to achieve weak recovery.
A.4 PROOF OF LEMMA 7
Denote Ã = aintran [ 1 2 (11 > + νν>)− I ] + ainter2n (11
>−νν>) as the expected adjacency matrix conditioned on community labels ν. Let Pij denote the set of all self-avoiding paths ik0 := {i0, . . . , ik} from i to j, such that no nodes appear twice in the path, and define Q{m}ij := ∑ im0 ∈Pij ∏m t=1(A − Ã)it−1it . Let T m ij be the concatenation of self-avoiding paths i k−m 0 and i k k−m+1, and let R m ij denote Tmij \ Pij . Define W {k,m} ij := ∑ ik0∈Rmij ∏k−m t=1 (A − Ã)it−1itÃik−mik−m+1 ∏k t=k−m+2(A)it−1it .
Then, by (Massoulié, 2014, Theorem 2.2), we can have
A{k} = Q{k} + k∑ m=1 (Q{k−m}ÃA{m−1})− k∑ m=1 W {k,m}.
Hence,
Γ {r} θ = r∑ k=r/2 θkA {k}
= r∑ k=r/2 θkQ {k} + r∑ k=r/2 θk k∑ m=1 ( Q{k−m}ÃA{m−1} ) − r∑ k=r/2 k∑ m=1 θkW {k,m}
= r∑ k=r/2 θkQ {k} + r∑ m=1  r∑ k=max(r/2,m) θkQ {k−m}Ã A{m−1} − r∑ k=r/2 k∑ m=1 θkW {k,m}.
In particular, for any u that is orthogonal toA{m}1 andA{m}ν with norm-1 (i.e., a feasible vector of the variational definition in Lemma 7), we have:∥∥∥Γ{r}θ u∥∥∥
2 ≤ ρ  r∑ k=r/2 θkQ {k} + r∑ m=1 r∑ k=max(r/2,m) |θk|ρ ( Q{k−m} )∥∥∥ÃA{m−1}u∥∥∥ 2
+ ρ  r∑ k=r/2 θk k∑ m=1 W {k,m}  Since by (Massoulié, 2014), the terms ρ ( Q{k} ) and ρ ( W {k,m}
) are bounded by n ξr/21 ,∥∥∥ÃA{m−1}x∥∥∥
2 is bounded by O(log(n) +
√ log(n)ξm−11 ), with some elementary calculations,
we have that ∥∥∥Γ{r}θ u∥∥∥
2 is bounded by ‖θ‖1n ξr/21 O(log(n)).
A.5 MORE EXPERIMENTAL DETAILS AND RESULTS
Experimental setup. We followed the setup of (Yang et al., 2016; Kipf & Welling, 2017) for citation networks Citeseer, Cora and Pubmed. The statistics of datasets is in Table A.1. Degree
distribution in the original graph and the powered graphs for Citeseer and Cora is shown in Figure A.1. We used the same dataset splits as in (Yang et al., 2016; Kipf & Welling, 2017) with 20 labels per class for training, 500 labeled samples for hyper-parameter setting, and 1, 000 samples for testing. To facilitate comparison, we use the same hyper-parameters as Vanilla GCN (Kipf & Welling, 2017), including number of layers (2), number of hidden units (16), dropout rates (0.5), weight decay (5× 10−4), and weight initialization following (Glorot & Bengio, 2010) implemented in Tensorflow. The training process is terminated if validation accuracy does not improve for 40 consecutive steps. We note that the citation datasets are extremely sensitive to initializations; as such, we report the test accuracy for the top 50 runs out of 100 random initializations sorted by the validation accuracy. We used Adam (Kingma & Ba, 2014) with a learning rate of 0.01 except for θ (VPN), which was chosen to be 1 × 10−5 to stabilize training. For VPN, we consider two-step training process, i.e., we first employ graph sparsification on adjacency matrix Ā[k] based on X for the first-step training, then process second graph sparsification on adjacency matrix Ā[k] with respect to the feature embedding from hidden layer and fetch them back into model for second-step training with warm restart (continue training with the parameters from first-step). We set αk to be 0.5 for the corresponding power order and 0 other wise in r-GCN.
To enhance robustness against evasion attack, we adaptively chose the threshold of graph sparsification. For each node, we consider all its neighbors within r hops, where r is the order of the powered graph. We order these neighbors based on their Euclidean distance to the current node in the feature space. We then set the threshold for this node such that only the first s∗d neighbors with the smallest distance are selected, where s is the sparsification rate, and d is the degree of the node in the original graph. If this node is a high-degree node (determined by if the degree is higher than the average degree by more than 2 standard deviations), we set s = 1; otherwise, we set s to be a small number slightly larger than 1. Nevertheless, we do not observe that the robustness performance is sensitive to this parameter. In addition, one can define their own distance function that better measures proximity between nodes, such as correlation distance or cosine distance. As a side remark, one can easily obtain A(r) with 1[(I + A)r > 0], and recursively calculate A[k]. The proposed methods have similar computational complexity as vanilla GCN as a result of sparsification.
Revisiting SBM. To obtain networks with different SNRs (0.58, 0.68, 0.79, 0.85, 0.91), we set the intra-degree average aintra = 2.1 and inter-degree average ainter to be 0.4, 0.3, 0.2, 0.15, 0.1 respectively. We generate 10 networks for each SNR. One of the networks with SNR of 0.91 is analyzed in Figure A.3 for different convolution operators. As predicted by Theorem 3, the proposed operators have nontrivial spectral gap. This is not enjoyed by other operators, so it is a quite unique property. Furthermore, our learned operator w/ or w/o normalization obtained 81% accuracy based
on the second eigenvector, whereas adjacency matrix catches “high-degree nodes” and normalized (or powered) Laplacians catch “tails”, with spectral clustering accuracy of only 51%.
Visualization of network embeddings. We used t-SNE (Maaten & Hinton, 2008) to visualize the network embeddings learned by r-GCN (Figure A.4) and VPN (Figure A.5) for Citeseer, Cora and Pubmed. As can be seen, with very limited labels, the proposed methods can learn an effective embedding of the nodes.
Evasion experiment. Five representative global adversarial attack methods are considered as comparisons:
• DICE (Zügner & Günnemann, 2019) (remove internally, insert externally): it randomly disconnects nodes from the same class and connects nodes from different classes.
• Aabr and ADW3 (Bojchevski & Günnemann, 2019): both attack methods are matrix perturbation theory based. ADW3 is the closed-form attack for edge deletion and Aabr is the add-by-remove strategy based method for edge insertion. • Meta-Train and Meta-Self (Zügner & Günnemann, 2019): both attack approaches are meta
learning based. The meta-gradient approach with self-training is denoted as Meta-Self and the variant without self-training is denoted as Meta-Train. Note that both attack methods result in OOM (out of memory) on Pubmed, we choose to report the performance against them only on Citeseer and Cora.
We provide the post-evasion accuracy in Table A.2 and Table A.3 and robustness enhancement in Figure A.6, Figure A.7 and Figure A.8 for Citeseer, Cora, and Pubmed. Results are obtained with 20 random initializations on a given attacked network.
Computational efficiency. The adjacency matrix could be much denser after being powered several times, which poses a general challenge for all high-order matrix-based approaches. To alleviate this issue, we employ a simple sparsification strategy as described in the main paper. To demonstrate the computational efficiency of the proposed method, we introduce a running time comparison experiment on the real world social network (Social circles: Facebook) (Leskovec & Mcauley, 2012). This small dataset consists of “circles” (or “friends lists”) from Facebook and becomes very dense from first order (4.71%) to second order (92.12%), thus is suitable to be used for computational efficiency evaluation. The results are shown in Table A.4. We can find that the running efficiency of the proposed method is compatible with baselines and the density does not affect the efficiency significantly when the dataset is not too large. Moreover, in this paper, our primary goal is to improve the robustness. We leave the solution to resolve the scalability issue as a future direction.

Review summaries:
This paper proposes a new architecture for graph convolutional network based on graph powering operation which generates a new graph based on the shortest distance between pair of nodes.  Its main motivation is to overcome the dominance of the first eigenvector in the existing GCN architectures based on the graph Laplacian operator. The theoretical evidence for the robustness is provided based on the signal-to-noise (SNR) ratio of the simplified stochastic block model (SBM). Two versions of the algorithms are proposed, namely the robust graph convolutional network (r-GCN) and variable power network (VPN). First, r-GCN is based on augmenting the graphs with graph powering operation. Next, VPN replaces the adjacency matrix of the graph convolutional operator by the newly proposed variable power operator. An additional sparsification scheme is proposed since the graph powering operation densifies the original graph.  

This paper proposes a graph convolutional operator based on graph powering and applies it to GCN architecture to improve the performance and robustness. This work is mainly motivated by the paper (Graph powering and spectral robustness, Abbe et al., 2018). The authors introduce the graph powering to graph convolution neural network domain to replace the original Laplacian operator. They further propose a graph sparsification/pruning strategy on the powered adjacency matrices in order to reduce the complexity and increase the robustness against adversarial attacks. They also provide theoretical analysis to prove that the proposed powering operator and subsequent methods have some spectral properties and theoretical feasibility. However, some conclusions are limited to the ideal situations or seem subjective. Extensive experiments are conducted to show better or comparable performance in both benign and adversarial situations.