Year: 2021

Abstract:
Empirical studies on the robustness of graph neural networks (GNNs) have suggested a relation between the vulnerabilities of GNNs to adversarial attacks and the increased presence of heterophily in perturbed graphs (where edges tend to connect nodes with dissimilar features and labels). In this work, we formalize the relation between heterophily and robustness, bridging two topics previously investigated by separate lines of research. We theoretically and empirically show that for graphs exhibiting homophily (low heterophily), impactful structural attacks always lead to increased levels of heterophily, while for graph with heterophily the change in the homophily level depends on the node degrees. By leveraging these insights, we deduce that a design principle identified to significantly improve predictive performance under heterophily—separate aggregators for egoand neighbor-embeddings—can also inherently offer increased robustness to GNNs. Our extensive empirical analysis shows that GNNs adopting this design alone can achieve significantly improved empirical and certifiable robustness compared to the best-performing unvaccinated model. Furthermore, models with this design can be readily combined with explicit defense mechanisms to yield improved robustness with up to 18.33% increase in performance under attacks compared to the best-performing vaccinated model.

heading: 1 INTRODUCTION
text: Graph neural networks (GNNs) aim to translate the enormous empirical success of deep learning to data defined on non-Euclidean domains, such as manifolds or graphs (Bronstein et al., 2017), and have become important tools to solve a variety of learning problems for graph structured and geometrically embedded data. However, recent works show that GNNs—much like their “standard” deep learning counterparts—have a high sensitivity to adversarial attacks: intentionally introduced minor changes in the graph structure can lead to significant changes in performance. This finding, first articulated by Zügner et al. (2018) and Dai et al. (2018), has triggered studies that investigated different attack scenarios (Xu et al., 2019; Wu et al., 2019; Li et al., 2020a; Ma et al., 2020).
A different aspect of GNNs that has been scrutinized recently is that most GNNs do not perform well with many heterophilous datasets. GNNs generally perform well under homophily (or assortativity), i.e., the tendency of nodes with similar features or class labels to connect (Pei et al., 2020; Zhu et al., 2020). Such datasets are thus called homophilous (or assortative). While homophilous datasets dominate the study of networks, homophily is not a universal principle; certain networks, such as romantic relationship networks or predator-prey networks in ecology, are mostly heterophilous (or disassortative). Employing a GNN which does not account for heterophily can lead to significant performance loss in heterophilous settings (Abu-El-Haija et al., 2019; Zhu et al., 2020; Bo et al., 2021). Previous works have thus proposed architectures for heterophilous data.
While previous work has focused on naturally-occurring heterophily, heterophilous interactions may also be introduced as adversarial noise: as many GNNs exploit homophilous correlation, they can be sensitive to changes that render the data more heterophilous. A natural follow-up question is if and how this observation manifests itself in previously proposed attacking strategies on GNNs. In this work, we thus investigate the relation between heterophily and robustness of GNNs against adversarial perturbations of graph structure, focusing on semi-supervised node classification tasks.
More specifically, our main contributions are:
• Formalization: We formalize the relation between adversarial structural attacks and the change of homophily level in the underlying graphs with theoretical (§3.1) and empirical (§5.1) analysis. Specifically, we show that on homophilous graphs, effective structural attacks lead to increased heterophily, while, on heterophilous graphs, they alter the homophily level contingent on node degrees. To our knowledge, this is the first formal analysis of such kind. • Heterophily-inspired Design: We show how the relation between attacks and heterophily can inspire more robust GNNs by demonstrating that a key architectural feature in handling heterophily, separate aggregators for ego- and neighbor-embeddings, also improves the robustness of GNNs against attacks (§3.2). • Extensive Empirical Analysis: We show the effectiveness of the heterophilous design in improving empirical (§5.2) and certifiable (§5.3) robustness of GNNs with extensive experiments on real-world homophilous and heterophilous datasets. Specifically, we compare GNNs with this design, which we refer to as heterophily-adjusted GNNs, to non-adjusted models, including stateof-the-art models designed with robustness in mind. We find that heterophily-adjusted GNNs are up to 5 times more certifiably robust and have stronger performance under attacks by up to 32.92% compared to non-adjusted, standard models. Moreover, this design can be combined with existing vaccination mechanisms, yielding up to 18.33% higher accuracy under attacks than the best non-adjusted vaccinated model.

heading: 2 NOTATION AND PRELIMINARIES
text: Let G = (V, E ,X) be a simple graph with node set V , edge set E and node attributes X. The one-hop neighborhood N(v) = {u : (u, v) ∈ E} of a node v ∈ V is the set of all nodes directly adjacent to v; the k-hop neighborhood of v ∈ V is the set of nodes reachable by a shortest path of length k. We represent the graph G algebraically by an adjacency matrix A ∈ {0, 1}|V|×|V| and node feature matrix X ∈ R|V|×F . We use As = A+ I to denote the adjacency matrix with self-loops added, and denote the corresponding row-stochastic matrices as Ā = D−1A and Ās = D−1s As, respectively, where D is a diagonal matrix with Dii = ∑ jAij (Ds is defined analogously). We further assume that there exists a vector y, which contains a unique class label yv for each node v. Given a training set TV = {(v1, y1), (v2, y2), ...} of labeled nodes, the goal of semi-supervised node classification is to learn a mapping ` : V → Y from the nodes to the set Y of class labels. Graph neural networks (GNNs). Most current GNNs operate according to a message passing paradigm where a representation vector rv is assigned to each node v ∈ V and continually updated by K layers of learnable transformations. These layers first aggregate representations over neighboring nodes N(v) and then update the current representation via an encoder ENC. For prevailing GNN models like GCN (Kipf & Welling, 2017) and GAT (Veličković et al., 2018), each layer can be formalized as r(k)v = ENC ( AGGR ({ r (k−1) u : u ∈ N(v) ∪ {v} })) , where AGGR is the mean func-
tion weighted by node degrees (GCN) or an attention mechanism (GAT), and ENC is a learnable (nonlinear) mapping.
Adversarial attacks on graphs. Given a graph G = (V, E ,X) and a GNN f that processes G, an adversarial attacker tries to create a perturbed graph G′ = (V, E ′,X) with a modified edge-set E ′ such that the performance of the GNN f is maximally degraded. The information available to the attacker can vary under different scenarios (Jin et al., 2020a; Sun et al., 2020). Here, we follow the gray-box formalization by Zügner et al. (2018), where the attacker knows the training set TV , but not the trained GNN f . The attacker thus considers a surrogate GNN and picks perturbations that maximize an attack loss Latk (Zügner et al., 2018; Zügner & Günnemann, 2019a), assuming that attacks to the surrogate model are transferable to the attacked GNN. For node classification, the attack loss Latk quantifies how the predictions zv ∈ [0, 1]|Y| made by the GNN f differ from the true labels y. For a targeted attack of node v with class label yv ∈ Y , we adopt the negative classification margin (CM-type) (Zügner et al., 2018; Xu et al., 2019): Latk = −∆c = −(zv,yv − maxy 6=yv zv,y). The attacker usually has additional constraints, such as a limit on the size of the perturbations allowed (Zügner et al., 2018; Zügner & Günnemann, 2019a).
Taxonomy of attacks. We follow the taxonomy of attacks introduced in (Jin et al., 2020a; Sun et al., 2020). For node classification, the attacker may aim to change the classification of a specific node v ∈ V (targeted attack), or to decrease the overall classification accuracy (untargeted attack). Attacks can also happen at different stages of the training process: we refer to attacks introduced
before training as (pre-training) poison attacks, and attacks introduced after the training process (and before potential retraining on perturbed data) as (post-training) evasion attacks. While our theoretical analysis (§3) mainly considers targeted evasion attacks, we consider other attacks in our empirical evaluation (§5).
Characterizing homophily and heterophily in graphs. Using class labels, we characterize the types of connections in a graph contributing to its overall level of homophily/heterophily as follows: Definition 1 (Homo/Heterophilous path and edge) A k-hop homophilous path from node w to u is a length-k path between endpoint nodes with the same class label yw = yu. Otherwise, the path is called heterophilous. A homophilous or heterophilous edge is a special case with k = 1.
Following (Zhu et al., 2020; Lim et al., 2021), we define the homophily ratio h as: Definition 2 (Homophily ratio) The homophily ratio is the fraction of homophilous edges among all the edges in a graph: h = |{(u, v) ∈ E|yu = yv}|/|E|. When the edges in a graph are wired randomly, independent to the node labels, the expectation for h is hr = 1/|Y| for balanced classes (Lim et al., 2021). For simplicity, we informally refer to graphs with homophily ratio h 1/|Y| as homophilous graphs (which have been the focus in most prior works), graphs with homophily ratio h 1/|Y| as heterophilous graphs, and graphs with homophily ratio h ≈ 1/|Y| as weakly heterophilous graphs.

heading: 3 RELATION BETWEEN GRAPH HETEROPHILY & MODEL ROBUSTNESS
text: In this section, we first show theoretical results on the relation between adversarial structural attacks and the change in the homophily level of the underlying graphs. Though empirical analyses from previous works have suggested this relation on homophilous graphs (Wu et al., 2019; Jin et al., 2020a), to our knowledge, we are the first to formalize it with theoretical analysis and address the case of heterophilous graphs. As an implication of the relation, we then discuss how a key design that improves predictive performance of GNNs under heterophily can also help boost their robustness.

heading: 3.1 HOW DO STRUCTURAL ATTACKS CHANGE HOMOPHILY IN GRAPHS?
text: Homophilous Graphs: Structural Attacks are Mostly Heterophilous Attacks. Our first result shows that, for homophilous data, effective structural attacks on GNNs (as measured by loss Latk) always result in a reduced level of homophily where either new heterophilous connections are added or existing homophilous connections are removed. It also states that direct perturbations on 1-hop neighbors of the target nodes are more effective than indirect perturbations (influencer attacks (Zügner et al., 2018)) on multi-hop neighbors. For simplicity, akin to previous works (Zügner et al., 2018; Zügner & Günnemann, 2019a) we establish our results for targeted evasion (post-training) attacks in a stylized learning setup with a linear GNN. However, our findings generalize to more general setups on real-world datasets as we show in our experiments (§5.1). In the theorems below, we use the notion of gambit node: node u is called a gambit if a perturbation that targets node v ∈ V adjusts the connectivity of node u ∈ V .
Theorem 1 Let G = (V, E ,X) be a self-loop-free graph with adjacency matrix A and node features xv = p · onehot(yv) + 1−p|Y| · 1 for each node v, where 1 is an all-1 vector, and p is a parameter that regulates the signal to noise ratio. Assume that a fraction h of each node’s neighbors belong to the same class, while a fraction 1−h|Y|−1 belongs uniformly to any other class. Consider a 2-layer linear GNN f (2)s (A,X) = Ā2sXW trained on a training set TV ⊆ DV , with at least one node from each class y ∈ Y , and degree d for all nodes with a distance less than 2 to any v ∈ DV . For a unit structural perturbation that involves a target node v ∈ DV , and a correctly classified gambit node with degree da, the following statements hold if h ≥ 1|Y| : 1. the attack loss Latk (§2) of the target v increases only for actions increasing heterophily, i.e.,
when removing a homophilous edge or path, or adding a heterophilous edge or path to node v; 2. direct perturbations on edges (or 1-hop paths) incident to the target node v lead to greater
increase in Latk than indirect perturbations on multi-hop paths to target node v.
We give the proof in App. C.1. Intuitively, the relative inability of existing GNNs to make full use of heterophilous data (Pei et al., 2020; Zhu et al., 2020) can be exploited by inserting heterophilous
connections in graphs where homophilous ones are expected. Though the theorem shows that effective attacks on homophilous graphs necessarily reduce the homophily level, the converse is not true: not all perturbations which reduce the homophily level are effective attacks (Ma et al., 2021).
Heterophilous Graphs: Structural Attacks Can Be Homophilous or Heterophilous, depending on Node Degrees. When a graph displays heterophily, our analysis shows a more complicated picture on how the level of homophily in the graph is changed by effective structural attacks: in heterophilous case, the direction of change is dependent on the degrees of both the target node v and the gambit node u of the attack. Specifically, if the degree of either node is low, attacks increasing the heterophily are still effective; however, if the degrees d and da of both nodes are high, attacks decreasing the heterophily will be effective. Similar to the homophilous case, we formalize our results below for targeted evasion attacks in a stylized learning setup. Theorem 2 Under the setup of Thm. 1, for a unit perturbation that involves a target node v with degree d, and a correctly classified gambit node with degree da, the following statements hold: 1. (Low-degree target node) if 0 < d ≤ |Y| − 2, for any da ≥ 0 and h ∈ [0, 1], the attack loss Latk (§2) of v increases only under actions increasing heterophily in the graph; 2. (High-degree target node) if d > |Y| − 2, conditioning on the degree da of the gambit node: (a) (Low-degree gambit node) if da < (d+2)(|Y|−1)d−|Y|+2 , for any h ∈ [0, 1], the attack loss Latk
(§2) of v increases only under actions increasing heterophily in the graph; (b) (High-degree gambit node) if da ≥ (d+2)(|Y|−1)d−|Y|+2 , for 0 ≤ h < da(d−|Y|+2)−(d+2)(|Y|−1) (d+1)|Y|da <
1 |Y| , Latk (§2) of v increases only under actions reducing heterophily in the graph.
In the statements above, the actions increasing heterophily in the graph include removing a homophilous edge or adding a heterophilous edge to node v, and the actions reducing heterophily in the graph include when adding a homophilous edge or removing a heterophilous edge to node v. The above theorems cover the situation when the gambit nodes are initially classified correctly (where attacks introducing heterophily can be unambiguously defined using the ground-truth class labels of the nodes involved). However, in §5.1, we show on real-world datasets that a relaxed interpretation of the theorems, where heterophily is instead defined by the predicted class labels of GNNs, can explain the behavior of the attacks regardless of the initial correctness of the gambits.

heading: 3.2 BOOSTING ROBUSTNESS WITH A SIMPLE HETEROPHILOUS DESIGN
text: A natural follow-up question is whether GNNs with better performance under heterophily are also more robust against structural attacks. We deduce that a key design for improving GNN performance for heterophilous data—separate aggregators for ego- and neighbor-embeddings—can also boost the robustness of GNNs by enabling them to better cope with adversarially-introduced changes in heterophily.
Separate Aggregators for Ego- and Neighbor-embeddings. This design uses separate GNN aggregators for ego-embedding rv and neighbor-embeddings {ru : u ∈ N(v)}. Formally, the representation learned for node v in the k-th layer is:
r(k)v = ENC ( AGGR1(r(k−1)v , r (k−2) v , ..., r (0) v ), AGGR2({r(k−1)u : u ∈ N(v)}) ) , (1)
where AGGR1 and AGGR2 are separate aggregators, such as averaging functions (GCN), attention mechanisms (GAT), or other pooling mechanisms (Hamilton et al., 2017). The ego-aggregator AGGR1 may also introduce skip connections (Xu et al., 2018) to the ego-embeddings aggregated in previous layers as shown in Eq. (1). This design has been utilized in existing GNN models (see §D.1 for details), and has been shown to significantly boost the representation power of GNNs under natural heterophily (Zhu et al., 2020).
Intuition. The key design changes, as compared to the GCN formulation in §2, allow for the egoembedding rv to be aggregated and weighted separately from the neighbor-embeddings {ru : u ∈ N(v)}, as well as for the use of skip connections to ego-embeddings of previous layers. Intuitively, ego-embeddings of feature vectors at the first layer are independent of the graph structure and thus unaffected by adversarial structural perturbations. Hence, a separate aggregator and skip connections can provide better access to unperturbed information, and helps mitigate the effects of the attacks.
Theoretical Analysis. We formalize the above intuition that shows how separate aggregators for ego- and neighbor-embeddings enable GNN layers to reduce the attack loss.
Theorem 3 Under the setup of Thm. 1, consider two alternative layers from which a two-layer linear GNN is built: (1) a layer defined as fs(A,X) = ĀsXW; and (2) a layer formulated as f(A,X;α) = ( (1− α)Ā + αI ) XW, which mixes the ego- and neighbor-embedding linearly under a predefined weight α ∈ [0, 1]. Then, for h > 1/|Y|, α > 1/(1 + da), and a unit perturbation increasing Latk as in Thm. 1, outputs of layer f lead to a strictly smaller increase in Latk than fs.
We provide the proof in App. C.3; note that for α = 1/(1 + da), the two layers are the same: f(A,X;α) = fs(A,X). Theorem 3 shows that an increase to the weights of ego-embedding improves the robustness of the GNN f for a homophily ratio h > 1/|Y|. Though aggregators and encoders are stylized in the theorem, the empirical analysis in §5.2 confirms that GNNs with more advanced aggregators and encoders also benefit from separate aggregators. Specifically, we find that such GNNs outperform methods without this design by up to 33.33% and 48.88% on homophilous and heterophilous graphs, respectively, while performing comparably on clean datasets.

heading: 4 RELATED WORK
text: Adversarial Attacks and Defense Strategies for Graphs Since NETTACK (Zügner et al., 2018) and RL-S2V (Dai et al., 2018) first demonstrated the vulnerabilities of GNNs against adversarial perturbations, a variety of attack strategies under different scenarios have been proposed, including adversarial attacks on the graph structure (Dai et al., 2018; Xu et al., 2019; Bojchevski & Günnemann, 2019a; Li et al., 2020a; Chang et al., 2020), node features (Takahashi, 2019; Ma et al., 2020), or combinations of both (Zügner et al., 2018; Zügner & Günnemann, 2019a; Wu et al., 2019). On the defense side, various techniques for improving the GNN robustness against adversarial attacks have been proposed, including: adversarial training (Xu et al., 2019; Zügner & Günnemann, 2019a; Bojchevski & Günnemann, 2019b); RGCN (Zhu et al., 2019), which adopts Gaussian-based embeddings and a variance-based attention mechanism; low-rank approximation of graph adjacency (Entezari et al., 2020) against Nettack (Zügner et al., 2018); Pro-GNN (Jin et al., 2020b), which estimates the unperturbed graph structure in training with the assumptions of low-rank, sparsity, and homophily of node features; GCN-Jaccard (Wu et al., 2019) and GNNGuard (Zhang & Zitnik, 2020), which assume homophily of features (or structural embeddings) and train GNN models on a pruned graph with only strong homophilous links; and Soft Medoid (Geisler et al., 2020), an aggregation function with improved robustness. Other recent works have looked into the certification of nodes that are guaranteed to be robust against certain structural and feature perturbations (Zügner & Günnemann, 2019b; Bojchevski & Günnemann, 2019b; Zügner & Günnemann, 2020), including approaches based on model-agnostic randomized smoothing (Cohen et al., 2019; Lee et al., 2019; Bojchevski et al., 2020). Interested readers are referred to the recent surveys (Jin et al., 2020a; Sun et al., 2020) for a comprehensive review of the literature.
GNNs & Heterophily Recent works (Pei et al., 2020; Liu et al., 2020; Zhu et al., 2020; Ma et al., 2021) have shown that heterophilous datasets can lead to significant performance loss for popular GNN architectures (e.g., GCN (Kipf & Welling, 2017), GAT (Veličković et al., 2018)). This issue is also known in classical semi-supervised learning (Peel, 2017). To address this issue, several GNN designs for handling heterophilous connections have been proposed (Abu-El-Haija et al., 2019; Pei et al., 2020; Zhu et al., 2020; Dong et al., 2021; Li et al., 2021; Zhu et al., 2021; Bo et al., 2021). Yan et al. (2021) recently discussed the connection between heterophily and oversmoothing for GNNs, and proposed designs to simultaneously address both issues. However, the formal connection between heterophily and robustness of GNNs has received little attention. Here we focus on a simple yet powerful design that significantly improves performance under heterophily (Zhu et al., 2020), and can be readily incorporated into GNNs.

heading: 5 EMPIRICAL EVALUATION
text: Our analysis seeks to answer the following questions: (Q1) Does our theoretical analysis on the relations between adversarial attacks and changes in heterophily level generalize to real-world datasets? (Q2) Do heterophily-adjusted GNNs, i.e., models with separate aggregators for ego- and neighborembeddings, show improved robustness against state-of-the-art attacks? (Q3) Does the identified design improve the certifiable robustness of GNNs?
First, we describe the experimental setup and datasets that we use to answer the above questions.
Attack Setup. We consider both targeted and untargeted attacks (§2), generated by NETTACK (Zügner et al., 2018) and Metattack (Zügner & Günnemann, 2019a), respectively. For each attack method, we consider poison (pre-training) and evasion (post-training) attacks, yielding 4 attack scenarios in total. We focus on robustness against structural perturbations and keep the node features unchanged. We randomly generate 3 sets of perturbations per attack method and dataset, and consistently evaluate each GNN model on them. We provide more details in App. D.2.
GNN Models. To show the effectiveness of our identified design, we evaluate four groups of models against adversarial attacks: (1) Models with this heterophilous design only: GraphSAGE (Hamilton et al., 2017), H2GCN (Zhu et al., 2020), CPGNN (Zhu et al., 2021), GPR-GNN (Chien et al., 2021) and FAGCN (Bo et al., 2021); we discuss how these models instantiate this design in App. §D.1; (2) State-of-the-art “vaccinated” architectures designed with robustness in mind: ProGNN (Jin et al., 2020b), GNNGuard (Zhang & Zitnik, 2020) and GCN-SVD (Entezari et al., 2020); (3) Models with both the heterophilous design and explicit robustness-enhancing mechanisms based on lowrank approximation: H2GCN-SVD and GraphSAGE-SVD; (4) Models without any vaccination, including some of the most popular methods: GCN (Kipf & Welling, 2017), GAT (Veličković et al., 2018), and the graph-agnostic multilayer perceptron (MLP) which relies only on node features. We discuss the implementations and parameters used for the models in App. D.
Datasets & Evaluation Setup. We consider two standard datasets with strong homophily, Cora (McCallum, 2000) and Citeseer (Sen et al., 2008), complemented with one weakly and one strongly heterophilous graph, introduced by Lim et al. (2021): FB100 (Traud et al., 2012) and Snap Patents (Leskovec et al., 2005; Leskovec & Krevl, 2014). We report summary statistics
in Table 1, and provide more details in App. D.5. For computational tractability, we subsample the Snap Patents data via snowball sampling (Goodman, 1961), where we keep 20% of the neighbors for each traversed node (see App. D.5 for details). We follow the evaluation procedure of Zügner et al. (2018) and Jin et al. (2020b), where we split the nodes of each dataset into training (10%), validation (10%) and test (80%) data, and determine the model parameters using the training and validation splits. We report the average performance and standard deviation on the 3 sets of generated perturbations. For targeted attacks with NETTACK, we report the classification accuracy on the target nodes; for untargeted attacks with Metattack, we report it over the whole test data.
Robustness Certificates. We adopt randomized smoothing for GNNs (Bojchevski et al., 2020) to evaluate the certifiable robustness, with parameter choices as detailed in App. D.2. We only consider structural perturbations in the randomization scheme. Following Geisler et al. (2020), we measure the certifiable robustness of GNN models with the accumulated certifications (AC) and the average maximum certifiable radii for edge additions (r̄a) and deletions (r̄d) over all correctly predicted nodes. More specifically, AC is defined as −R(0, 0) + ∑ ra,rd≥0R(ra, rd), where R(ra, rd) is the certifiably correct ratio, i.e., the ratio of the nodes in the test splits that are both predicted correctly by the smoothed classifier and certifiably robust at radius (ra, rd). In addition, we report the accuracy of each model with randomized smoothing enabled on the test splits of the clean datasets, which is equal to R(0, 0). We report the average and standard deviation of each statistic over the 3 different training, validation and test splits.

heading: 5.1 (Q1) STRUCTURAL ATTACKS ARE MOSTLY HETEROPHILOUS: EMPIRICAL VALIDATION
text: To show that our theoretical analysis in §3.1 generalizes to more complex settings beyond the assumptions we made in the theorems, we look into effective targeted attacks made by NETTACK on real-world homophilous and heterophilous datasets, and present statistics of the attacks in Table 2, with a focus on the ratios of heterophilous attacks. We use a budget of 1 perturbation per target node in this experiment, and the statistics are reported among all effective perturbations targeting nodes that are correctly classified on clean datasets by the surrogate GNN used by NETTACK (i.e., GCN) as described in §5. To help validate the dependency between the degrees of the target/gambit nodes
and the changes of heterophily predicted by Thm. 2, we also show the scatter plots of target and gambit node degrees in Fig. 2 in App. E.5.
Homophilous Networks. For the strongly homophilous Cora and Citeseer graphs, all changes in the graph structure that are introduced by effectives attacks follow the conclusion of Thm. 1: they reduce homophily (increase heterophily) by adding heterophilous edges or removing homophilous edges. These results show that despite the simplified analysis, the takeaway of Thm. 1 can be generalized to real-world datasets. In addition, the attacks mostly introduce, rather than prune, edges, suggesting that attacks adding outlier edges to the graph are more powerful
than attacks removing informative existing edges. These observations in our experiments are consistent with the observations from previous works (Jin et al., 2020a; Geisler et al., 2020).
Heterophilous Networks. For heterophilous graphs FB100 (h ≈ 1/|Y|) and Snap (h < 1/|Y|), Fig. 2 in App. E.5 shows that almost all attacks leverage gambit nodes with low degrees (1 or 2); no node with degree higher than 5 is leveraged. All attacks leveraging correctly classified gambit nodes are connecting node u ∈ V with a different ground-truth class label yu 6= yv to the target nodes v ∈ V; attacks leveraging incorrectly classified gambit nodes are always connecting node u with a different predicted class label ŷu 6= ŷv = yv to the target node v, even though some gambit nodes have the same ground-truth class label yu = yv 6= ŷu as the target nodes. These results validate the conclusion of Thm. 2 on correctly classified gambit nodes, and demonstrate its generalizability under the heterophily definition based on predicted class labels. Note that the predicted class labels ŷu for each node u ∈ V are based on GCN, which is the surrogate GNN used by NETTACK.

heading: 5.2 (Q2) BENCHMARK STUDY OF GNN MODELS: HETEROPHILOUS DESIGN LEADS TO IMPROVED EMPIRICAL ROBUSTNESS
text: To answer (Q2) on whether heterophily-adjusted GNN models show improved performance against state-of-the-art attacks, we conduct a comprehensive benchmark study. We consider all four categories of GNN models mentioned in §5, and evaluate their robustness against both targeted and untargeted attacks. We report the hyperparameters for each method in App. D.4. Table 3 shows the performance of each method under poison (pre-training) attacks, and Fig. 1 visualizes the corresponding performance changes relative to the clean datasets. For conciseness, we report further results under evasion (post-training) attacks and on clean (unperturbed) data in the Appendix (Table 7 for NETTACK; Table 8 for Metattack). Also, in App. E.4 we discuss how our simple heterophilous design leads to only minor computational overhead compared to existing vaccination mechanisms.
Targeted attacks by NETTACK. 1© Poison attacks. Under targeted poison attacks, Table 3 (left) shows that GraphSAGE-SVD and H2GCN-SVD, which combine our identified design with a lowrank vaccination approach adopted in GCN-SVD (Entezari et al., 2020) (details in App. D.3), outperform state-of-the-art methods across all datasets by up to 13.34% in homophilous settings and 18.33% in heterophilous settings.
Methods merely employing the identified design also show significantly improved robustness, though there are differences in the amount of robustness improvement due to architectural differences. Specifically, these methods outperform the best unvaccinated method (GAT) on all datasets by up to 32.92% in average, despite having comparable performance on clean datasets. On Citeseer and FB100, methods with the heterophilous design also show comparable or even better robustness than state-of-the-art vaccinated GNNs like ProGNN and GCN-SVD.
2© Evasion attacks. Under evasion attacks (Table 7), we observe similar trends as in poison attacks: GraphSAGE-SVD and H2GCN-SVD are up to 20.55% more accurate than the GCN-SVD baseline, and methods featuring the identified design alone achieve up to 24.45% gain in average performance against the best unvaccinated baseline. We note that GNNGuard and ProGNN are not capable of ad-
dressing evasion attacks. In summary, these observations show that the heterophily-inspired design is orthogonal to existing vaccination mechanisms for improving the robustness of GNNs.
Untargeted attacks by Metattack. 1© Poison attacks. We also test the robustness of each method against untargeted attacks. Table 3 (right) shows the performance under poison attacks. Though our theoretical analysis in §3 focuses on the effect of the heterophilous design under targeted attacks, we observe similar improvements in robustness against untargeted attacks in the poison setup. GNNs with the identified design show mostly improved robustness compared to unvaccinated models, while having similar performance on the clean datasets. Specifically, CPGNN shows exceptional robustness, outperforming the best unvaccinated model by up to 32.85%. Moreover, models combining
the identified design with low-rank approximation show more than 10% improvement accuracy compared to GCN-SVD, which uses only low-rank approximation, and ProGNN. We note that GNNGuard, which uses a similarity-based defense strategy, also shows more competitive robustness under untargeted attacks. This is in line with existing works showing that low-rank based approach does not adapt well to untargeted attacks (Jin et al., 2020b). Nevertheless, methods combining a heterophilous design and low-rank approximations still outperform GNNGuard on two datasets. 2© Evasion attacks. We present the performance under evasion attacks and on clean datasets in Table 8 in the Appendix. Unlike the poison attacks, the evasion setup only leads to a slight decrease in average accuracy of less than 2% for most models. Moreover, there appears to be no clearly increased robustness for vaccinated models (with the identified design or other vaccination machanisms) compared to unvaccinated models. This can be attributed to the reduced effectiveness of evasion vs. poison attacks (as in NETTACK), and the increased challenges of untargeted attacks.

heading: 5.3 (Q3) HETEROPHILY-ADJUSTED GNNS ARE CERTIFIABLY MORE ROBUST
text: It is worth noting that robustness against specific attacks such as NETTACK and Metattack does not guarantee robustness towards other possible attacks. To overcome this limitation, robustness certificates provide guarantees (in some cases probabilistically) that attacks within a certain radius cannot change a model’s predictions. Complementary to our evaluation on empirical robustness, we further demonstrate that heterophily-adjusted GNNs featuring our identified design are certifiably
more robust than methods without it, thus answering (Q3). For GNN models, we include H2GCN, GraphSAGE, CPGNN, GPR-GNN, FAGCN, GAT and GCN in this analysis. We exclude other models that either learn to rewrite the graph structure through the training process or require a recalculation of the low-rank approximation for every randomized perturbation as sampling on these models is infeasible. We use the same hyperparameters as in the benchmark study in §5.2.
Table 4 shows multiple metrics of certifiable robustness of each GNN model under a randomization scheme allowing for both addition and deletion of edges; we additionally report results under randomization schemes allowing only addition or deletion in Table 9. For the scheme allowing both addition and deletion, we observe that all heterophily-adjusted methods have better certifiable robustness compared to methods without the design. Specifically, on homophilous datasets (Cora and Citeseer), methods with the identified design achieve an up to 1.6 times relative improvement in accumulated certification. On heterophilous datasets (FB100 and Snap), they outperform the baselines by a factor of 4.4. In the more challenging case with the addition only scheme, methods with the design also show up to 1.2 times relative increase in AC on the homophilous datasets and 5.0 times relative increase in AC on the heterophilous datasets compared to the baselines. For the deletion only scheme, unvaccinated models like GCN already have decent certifiable robustness in this scenario. This is commensurate with our discussions in §5.1 that deletions create less severe perturbations. Overall, our results show that models featuring our identified design achieve significantly improved certifiable robustness compared to models lacking this design. However, like in our empirical robustness evaluation, architectural differences lead to some variability of robustness; the results also show tradeoffs between accuracy and robustness. We also observe that the robustness rankings for methods under certifiable and empirical robustness are different, as in the previous results from Geisler et al. (2020); we discuss possible reasons in App. E.3.

heading: 6 CONCLUSION
text: We formalized the relation between heterophily and adversarial structural attacks, and showed theoretically and empirically that effective attacks gravitate towards increasing heterophily in both homophilous and heterophilous graphs by leveraging low-degree (gambit) nodes. Using these insights, we showed that a key design addressing heterophily, namely separate aggregators for ego- and neighbor-embeddings, can lead to competitive improvement on empirical and certifiable robustness, with only small influence on clean performance. Finally, we compared the design with state-of-theart vaccination mechanisms under different attack scenarios for various datasets, and illustrated that they are complementary and that their combination can lead to more robust GNN models. We note that while we focus on the structural attacks, GNNs are also vulnerable to other types of attacks such as feature perturbations. In addition, the graph-agnostic MLP, which is immune to structural attacks, on some datasets outperforms all GNNs against attacks; this calls for further works in understanding of the nature of the attacks and also for effective defense strategies upon our discoveries.

heading: REPRODUCIBILITY STATEMENT
text: For reproducibility, we describe the experimental setups in §5 with additional details in App. D.2, including the implementation and the hyperparameter settings used for each method, and the hardware specifications. We provide anonymized code and datasets in the supplementary material with instructions for replicating the experiments, and will make them publicly available upon acceptance.

heading: A SOCIETAL IMPACTS
text: A large number of popular GNN models are inherently based on homophily. Our work shows that such models may be less robust to adversarial perturbations, and thus when employed for decision-making these models may lead to undesirable, erroneous or biased results. For example, an inherently homophilous GNN model may lead to the so-called “filter bubble” phenomenon in a recommendation system in which existing beliefs or preferences are reinforced. Similarly, as homophily-based GNNs typically average over node neighborhoods, this may result in less visibility of minority groups in the network, thus reinforcing disparities. Heterophilous network design may improve some of these aspects, and also benefit from some additional robustness as shown in this work.
However, it should be noted that heterophilous GNNs on their own cannot fully solve the aforementioned issues. In particular, while heterophilous design can improve robustness, it is not (in general) constructed to ensure other important aspects such as fairness. Better understanding of GNNs and tailored auditing tools are necessary in order to deploy these learning algorithms in the context of decision-making affecting humans.

heading: B NOMENCLATURE
text: We summarize the main symbols used in this work and their definitions below:

heading: C PROOFS AND DISCUSSIONS OF THEOREMS

heading: C.1 DETAILED ANALYSIS OF THEOREM 1
text: Proof 1 (for Thm. 1) We give the proof in three parts: first, we analyze the training process of the GNN f (2)s (A,X) = Ā2sXW on clean data and analytically derive the optimal weight matrix W∗ in a stylized learning setup; then, we construct a targeted evasion attack and calculate the attack loss for a unit structural attack; last, we summarize and validate the statements in the theorem.
Stylized learning on clean data. Given the 2-layer linearized GNN f (2)s (A,X) = Ā2sXW and the training set TV ⊆ DV , the goal of the training process is to optimize the weight matrix W to minimize the cross-entropy loss function L([z]TV ,:, [Y]TV ,:), where predictions [z]TV ,: = [Ā2sX]TV ,:W correspond to the predicted class label distributions for each node v in the training set TV , and [Y]TV ,: is the one-hot encoding of class labels provided in the training set.
Without loss of generality, we reorder TV accordingly such that the one-hot encoding of labels for nodes in the training set [Y]TV ,: is in increasing order of the class label yv:
[Y]TV ,: =  1 0 0 · · · 0 ... ... ... . . . ... 1 0 0 · · · 0 0 1 0 · · · 0 ... ... ... . . . ... 0 1 0 · · · 0 ... ... ... . . . ... 0 0 0 · · · 1 ... ... ... . . . ...
0 0 0 · · · 1  |TV |×|Y|
(2)
Now we look at the term [Ā2sX]TV ,: in [z]TV ,: = [Ā 2 sX]TV ,:W, which are the feature vectors aggregated by the two GNN layers for nodes v in the training set TV . As stated in the theorem, we assume TV ⊆ DV , where node u ∈ DV have degree d; proportion h of their neighbors belong to the same class, while proportion 1−h|Y|−1 of them belong to any other class uniformly, and for each node v ∈ V the node features are given as xv = p · onehot(yv) + 1−p|Y| · 1 for each node v ∈ V . Then, after the first layer, we have:
[ĀsX]TV ,: = 1
d+ 1  (hd+ 1)p 1−h|Y|−1dp 1−h |Y|−1dp · · · 1−h |Y|−1dp ... ... ... . . . ... (hd+ 1)p 1−h|Y|−1dp 1−h |Y|−1dp · · · 1−h |Y|−1dp 1−h |Y|−1dp (hd+ 1)p 1−h |Y|−1dp · · · 1−h |Y|−1dp ... ... ... . . . ... 1−h |Y|−1dp (hd+ 1)p 1−h |Y|−1dp · · · 1−h |Y|−1dp ... ... ... . . . ... 1−h |Y|−1dp 1−h |Y|−1dp 1−h |Y|−1dp · · · (hd+ 1)p ... ... ... . . .
... 1−h |Y|−1dp 1−h |Y|−1dp 1−h |Y|−1dp · · · (hd+ 1)p  |TV |×|Y| + 1− p |Y|
(3)
and after the second layer
[Ā2sX]TV ,: = 1
(d+ 1)2|Y|(|Y| − 1)  S1 T1 T1 · · · T1 ... ... ... . . . ... S1 T1 T1 · · · T1 T1 S1 T1 · · · T1 ... ... ... . . . ... T1 S1 T1 · · · T1 ... ... ... . . . ... T1 T1 T1 · · · S1 ... ... ... . . . ...
T1 T1 T1 · · · S1  |TV |×|Y| + 1 |Y| , (4)
where S1 = ((h|Y| − 1)d+ |Y| − 1)2 p, and T1 = ((h|Y|−1)d+|Y|−1) 2p
|Y|−1 .
For [Y]TV ,: and [Ā 2 sX]TV ,: which we derived in Eq. (2) and (4), we can find the optimal weight matrix W∗ such that [Ā2sX]TV ,:W∗ = [Y]TV ,:, making the cross-entropy loss L([z]TV ,:, [Y]TV ,:) = 0.
To find W∗, we can proceed as follows. First, sample one node from each class to form a smaller set TS ⊂ TV . Therefore, we have:
[Y]TS ,: =  1 0 0 · · · 0 0 1 0 · · · 0 ... ... ... . . . ...
0 0 0 · · · 1  |Y|×|Y| = I|Y|×|Y|
and
[Ā2sX]TS ,: = 1
(d+ 1)2|Y|(|Y| − 1)  S1 T1 T1 · · · T1 T1 S1 T1 · · · T1 ... ... ... . . . ...
T1 T1 T1 · · · S1  |Y|×|Y| + 1 |Y| .
Note that [Ā2sX]TS ,: is a specific circulant matrix, and therefore its inverse exists. Using the Sherman-Morrison formula, we can find its inverse as:
( [Ā2sX]TS ,: )−1 =
(d+ 1)2(|Y| − 1)2
p(d(h|Y| − 1) + |Y| − 1)2|Y| ·
 |Y| − 1 −1 · · · −1 −1 |Y| − 1 · · · −1
... ...
. . . ...
−1 −1 · · · |Y| − 1
+ 1
|Y|  . (5) Now, let W∗ = ( [Ā2sX]TS ,: )−1 , then we have [z]TS ,: = [Ā 2 sX]TS ,:W∗ = [Y]TS ,: = I|Y|×|Y|. It is also easy to verify that [z]TV ,: = [Ā 2 sX]TV ,:W∗ = [Y]TV ,:. Since W∗ satisfiesL([z]TV ,:, [Y]TV ,:) =
0, we know W∗ = ( [Ā2sX]TS ,: )−1 is the optimal weight matrix that we can learn under TV .
Attack loss under evasion attacks. Now consider an arbitrary target node v ∈ DV with class label yv ∈ Y , and a unit structural perturbation leveraging gambit node u ∈ V with degree da that affects the predictions zv of node v made by GNN f (2) s . Without loss of generality, we assume node v has yv = 1. As f (2) s contains 2 GNN layers with each layer aggregating feature vectors within neighborhoodN(v) of each node v, the perturbation must take place in the direct (1-hop) neighborhood N(v) or 2-hop neighborhood N2(v) to affect the predictions zv . For the unit perturbation, the
attacker can add or remove a homophilous edge or path between nodes u and v, which we denote as δ1 (δ1 = 1 for addition and δ1 = −1 for removal); alternatively, the attacker can add or remove a heterophilous edge or path between nodes u and v, which we denote as δ2 = ±1 analogously. We denote the perturbed graph adjacency matrix as Ā′s, and z ′ v = [Ā ′2 s X]v,:W∗
1© Unit perturbation in direct neighborhoodN(v). We first consider a unit perturbation in the direct (1-hop) neighborhoodN(v) of node v. For simplicity of derivation, we assume that the perturbation does not change the row-stochastic normalization of Ās, and only affects the aggregated feature vectors of the target node v.
In the case of δ1 = ±1 and δ2 = 0, we have
[Ā′sX]v,: − [ĀsX]v,: = 1
da + 1
[ δ1 ( 1−p |Y| + p ) δ1 ( 1−p |Y| ) · · · δ1 ( 1−p |Y| ) ] and
[Ā′2s X]v,: − [Ā2sX]v,: = 1
(da + 1)2(d+ 1)|Y|
[ S2 T2 |Y|−1 · · · T2 |Y|−1 ] ,
where S2 = δ1 (da(p(d(h|Y| − 1) + h|Y|+ |Y| − 2) + d+ 2) + (d+ 2)(|Y| − 1)p+ d+ 2) and T2 = δ1 (− (da(p(d(h|Y| − 1) + h|Y|+ |Y| − 2) + (−d− 2)(|Y| − 1)) + (d+ 2)(|Y| − 1)(p− 1))). By Multiplying [Ā′2s X]v,: by W∗, we can get the predictions z ′ v after perturbations; we omit the analytical expression of z′v here due to its complexity.
On the perturbed graph, the CM-type attack loss is calculated as
LCMatk (z′v) = −(z′v,yv −maxy 6=yv z′v,y)
Since on clean data, LCMatk (zv) = −1, the change in attack loss before and after attack is
∆LCMatk = LCMatk (z′v)− LCMatk (zv) = LCMatk (z′v) + 1 (6)
= − (d+ 1)δ1(|Y| − 1) (da(d(h|Y| − 1) + h|Y|+ |Y| − 2) + (d+ 2)(|Y| − 1)) (da + 1)2(d(h|Y| − 1) + |Y| − 1)2 . (7)
Solving following system of inequalities for δ1, ∆LCMatk > 0 h ∈ [0, 1] |Y| ≥ 2 d, da, |Y| ∈ Z+
(8)
we get the valid range of δ1 as δ1 < 0, when 0 < d ≤ |Y| − 2 δ1 < 0, when d > |Y| − 2 and da < (d+2)(|Y|−1)d−|Y|+2 δ1 < 0, when d > |Y| − 2 and da ≥ (d+2)(|Y|−1)d−|Y|+2 and 1 ≥ h > da(d−|Y|+2)−(d+2)(|Y|−1) (d+1)|Y|da
δ1 > 0, when d > |Y| − 2 and da ≥ (d+2)(|Y|−1)d−|Y|+2 and 0 ≤ h < da(d−|Y|+2)−(d+2)(|Y|−1) (d+1)|Y|da
.
(9) Note that the above solution is not applicable when h = d−|Y|+1d|Y| , in which case d(h|Y|−1)+ |Y|− 1 = 0 and the solution of optimal weight matrix W∗ = ( [Ā2sX]TS ,: )−1 is undefined.
In the case of δ1 = 0 and δ2 = ±1, we have
[Ā′sX]v,: − [ĀsX]v,: = 1
da + 1
[ δ2 ( 1−p |Y| ) δ2 ( 1−p |Y| + p ) · · · δ2 ( 1−p |Y| + p ) ] and
[Ā′2s X]v,: − [Ā2sX]v,: = 1
(da + 1)2(d+ 1)|Y|
[ S4 |Y|−1 T4 · · · T4 ] ,
where S4 = δ2 (− (da(p(d(h|Y| − 1) + h|Y|+ |Y| − 2) + (−d− 2)(|Y| − 1)) + (d+ 2)(|Y| − 1)(p− 1))) and T4 = δ2 (da(p(d(h|Y| − 1) + h|Y|+ |Y| − 2) + d+ 2) + (d+ 2)(|Y| − 1)p+ d+ 2). By
multiplying [Ā′2s X]v,: with W∗, we can get the predictions z ′ v after the perturbations. Following a similar derivation to that in the previous case, we can compute the change in the CM-type attack loss before and after attack as
∆LCMatk = (d+ 1)(|Y| − 1) (da(d(h|Y| − 1) + h|Y|+ |Y| − 2) + (d+ 2)(|Y| − 1)) δ2
(da + 1) 2(d(h|Y| − 1) + |Y| − 1)2 . (10)
Solving the same system of inequalities as Eq. (8) for δ2, we obtain the valid range of δ2 as δ2 > 0, when 0 < d ≤ |Y| − 2 δ2 > 0, when d > |Y| − 2 and da < (d+2)(|Y|−1)d−|Y|+2 δ2 > 0, when d > |Y| − 2 and da ≥ (d+2)(|Y|−1)d−|Y|+2 and 1 ≥ h > da(d−|Y|+2)−(d+2)(|Y|−1) (d+1)|Y|da
δ2 < 0, when d > |Y| − 2 and da ≥ (d+2)(|Y|−1)d−|Y|+2 and 0 ≤ h < da(d−|Y|+2)−(d+2)(|Y|−1) (d+1)|Y|da
.
(11) Note that the above solution is not applicable when h = d−|Y|+1d|Y| , in which case d(h|Y|−1)+ |Y|− 1 = 0 and the solution of optimal weight matrix W∗ = ( [Ā2sX]TS ,: )−1 is undefined. We also note that we always have da(d−|Y|+2)−(d+2)(|Y|−1)(d+1)|Y|da − 1 |Y| = − (|Y|−1)(da+d+2) (d+1)|Y|da < 0 when |Y| ≥ 2.
2© Unit perturbation in 2-hop neighborhood N2(v). We now consider a unit perturbation in the 2-hop neighborhood N(v) of node v. In this case we will have [Ā′sX]v,: = [ĀsX]v,:.
In the case of δ1 = ±1 and δ2 = 0, we have
[Ā′2s X]v,: − [Ā2sX]v,: = 1
(da + 1)2|Y|
[ S5 T5 |Y|−1 · · · T5 |Y|−1 ] ,
where S5 = (da(p(h|Y| − 1) + 1) + (|Y| − 1)p+ 1) δ1
and T5 = (da(−h|Y|p+ |Y|+ p− 1) + |Y|(−p) + |Y|+ p− 1) δ1. By multiplying [Ā′2s X]v,: with W∗, we can get the predictions z ′ v after perturbations. Following a similar derivation as before, we can get the change in the CM-type attack loss before and after attack as
∆LCMatk = − (d+ 1)2(|Y| − 1) (da(h|Y| − 1) + |Y| − 1) δ1
(da + 1) 2(d(h|Y| − 1) + |Y| − 1)2 . (12)
Solving the same system of inequalities as Eq. (8) for δ1, we get the valid range of δ1 as δ1 < 0, when da < |Y| − 1 δ1 < 0, when da ≥ |Y| − 1 and da−|Y|+1|Y|da < h ≤ 1 δ1 > 0, when da > |Y| − 1 and 0 ≤ h < da−|Y|+1|Y|da . (13)
Note that the above solution is not applicable when h = d−|Y|+1d|Y| , in which case d(h|Y|−1)+ |Y|− 1 = 0 and the solution of optimal weight matrix W∗ = ( [Ā2sX]TS ,: )−1 is undefined.
For the case δ1 = 0 and δ2 = ±1, we have
[Ā′2s X]v,: − [Ā2sX]v,: = 1
(da + 1)2|Y|
[ S6 |Y|−1 T6 · · · T6 ] where S6 = δ2 (da(−h|Y|p+ |Y|+ p− 1) + |Y|(−p) + |Y|+ p− 1) and T6 = δ2 (da(p(h|Y| − 1) + 1) + (|Y| − 1)p+ 1) . Multiplying [Ā′2s X]v,: with W∗, we can get the predictions z ′ v after perturbations. As before we can compute the change in the CM-type attack loss before and after attack as
∆LCMatk = (d+ 1)2(|Y| − 1) (da(h|Y| − 1) + |Y| − 1) δ2
(da + 1) 2(d(h|Y| − 1) + |Y| − 1)2 (14)
Finally, solving the same system of inequalities as Eq. (8) for δ2, we get the valid range of δ2 as δ2 > 0, when da < |Y| − 1 δ2 > 0, when da ≥ |Y| − 1 and da−|Y|+1|Y|da < h ≤ 1 δ2 < 0, when da > |Y| − 1 and 0 ≤ h < da−|Y|+1|Y|da . (15)
Note that the above solution is not applicable when h = d−|Y|+1d|Y| , in which case d(h|Y|−1)+ |Y|− 1 = 0 and the solution of optimal weight matrix W∗ = ( [Ā2sX]TS ,: )−1 is undefined. We also note that we always have da−|Y|+1|Y|da − 1 |Y| = 1−|Y| |Y|da < 0 when |Y| ≥ 2.
Summary and validation of theorem statements. Based on our derivations, we summarize and validate our statements in the theorem next.
1© The attack losses Latk (CM-type, §2) increase only by removing a homophilous edge or path, or adding a heterophilous edge or path to node v. From Eq. (9), (11), (13), (15), we observe that for both direct attacks in 1-hop neighborhood N(v) and indirect attacks in 2-hop neighborhood N2(v), when h ≥ 1|Y| , the attack loss increases only if δ1 < 0, which represents removal of a homophilous edge or path to node v, or if δ2 > 0, which represents addition of a heterophilous edge or path to node v.
2© Direct perturbations on edges (or 1-hop paths) of the target node v lead to greater increase in Latk than indirect perturbations on multi-hop paths to target node v. From Eq. (7) and Eq. (10), the change in the CM-type attack loss ∆LCMatk (z′v) for direct perturbations on 1-hop neighborhood N(v) of the target node v considering both δ1 and δ2 can be written as
∆LCM,directatk = (d+ 1)(|Y| − 1) (da(d(h|Y| − 1) + h|Y|+ |Y| − 2) + (d+ 2)(|Y| − 1)) (δ2 − δ1)
(da + 1) 2(d(h|Y| − 1) + |Y| − 1)2 (16)
From Eq. (12) and Eq. (14), the change in the CM-type attack loss ∆LCMatk (z′v) for indirect perturbations on 2-hop neighborhood N2(v) of the target node v considering both δ1 and δ2 is
∆LCM,indirectatk = (d+ 1)2(|Y| − 1) (da(h|Y| − 1) + |Y| − 1) (δ2 − δ1)
(da + 1) 2(d(h|Y| − 1) + |Y| − 1)2 (17)
Note that when h ≥ 1|Y| , we have h|Y| ≥ 1, and
dh|Y| − d+ 2|Y| − 2 d(h|Y| − 1) + |Y| − 1 = 1 + da(d(h|Y| − 1) + |Y| − 1) + (d+ 1)(|Y| − 1) da(h|Y| − 1) + |Y| − 1 > 1 (18)
Therefore we will always have ∆LCM,directatk > ∆L CM,indirect atk for an effective unit perturbation that increases attack loss LCMatk (i.e., δ1 = −1 and δ2 = 0, or δ1 = 0 and δ2 = 1) when h ≥ 1|Y| .

heading: C.2 DETAILED ANALYSIS OF THEOREM 2
text: Proof 2 (for Thm. 2) For a direct unit perturbation in the 1-hop neighborhood N(v) of the target node v, from Eq. (9) and (11) of Proof 1, we observe that the signs of δ1 and δ2 which increase the attack loss are contingent on the degree of the target node d, the degree of the gambit node da and the homophily ratio h of the graph:
1© if 0 < d ≤ |Y| − 2 (i.e., when degree d of the target node is low), regardless of da and h, the attack loss increases only if δ1 < 0, which represents removal of a homophilous edge to node v, or if δ2 > 0, which represents addition of a heterophilous edge to v;
2© if d > |Y| − 2 (i.e., when degree d of the target node is high), the increase of the attack loss will be dependent to the degree of the gambit node da and the homophily ratio h of the graph:
(a) when da < (d+2)(|Y|−1) d−|Y|+2 (i.e., when degree da of the gambit node is low), regardless of h, the
attack loss increases only if δ1 < 0, which represents removal of a homophilous edge to node v, or if δ2 > 0, which represents addition of a heterophilous edge to v;
(b) when da ≥ (d+2)(|Y|−1)d−|Y|+2 (i.e., when degree da of the gambit node is high), for 0 ≤ h < da(d−|Y|+2)−(d+2)(|Y|−1)
(d+1)|Y|da < 1 |Y| , the attack loss increases only if δ1 < 0, which represents
removal of a homophilous edge to node v, or if δ2 > 0, which represents addition of a heterophilous edge to v.
We note that the above conclusions are not applicable when h = d−|Y|+1d|Y| , in which case d(h|Y| − 1) + |Y| − 1 = 0 and the solution of optimal weight matrix W∗ = ( [Ā2sX]TS ,: )−1 is undefined.

heading: C.3 DETAILED ANALYSIS OF THEOREM 3
text: Proof 3 (for Thm. 3) In this proof, we mainly focus on analyzing the increase in Latk for the GNN layer defined as f(A,X;α) = ( (1− α)Ā + αI ) XW. We follow a similar process as in Proof 1,
since the layer defined as fs(A,X) = ĀsXW is a special case of the previous formulation when α = 11+d .
Layer f(A,X;α) = ( (1− α)Ā + αI ) XW. We first derive the optimal weight matrix W∗ in a stylized learning setup as in Proof 1. Following a similar process, for this GNN layer we have
[( (1− α)Ā + αI ) X ] TS ,: = 1− p |Y| +  S7 T7 T7 · · · T7 T7 S7 T7 · · · T7 ... ... ... . . . ...
T7 T7 T7 · · · S7  |Y|×|Y|
where S7 = p(α+ h− αh), and T7 = (α−1)(h−1)p|Y|−1 , and
[Y]TS ,: =  1 0 0 · · · 0 0 1 0 · · · 0 ... ... ... . . . ...
0 0 0 · · · 1  |Y|×|Y| = I|Y|×|Y|
Using the Sherman-Morrison formula, we find its inverse:([( (1− α)Ā + αI ) X ] TS ,: )−1 = |Y| − 1 p(a(h− 1)|Y| − h|Y|+ 1)|Y|
· 1− |Y| 1 1 · · · 1 1 1− |Y| 1 · · · 1 ... ... ... . . . ...
1 1 1 · · · 1− |Y| + 1|Y| Assuming W∗ = ([( (1− α)Ā + αI ) X ] TS ,: )−1 , we obtain
[z]TS ,: = [( (1− α)Ā + αI ) X ] TS ,: W∗ = [Y]TS ,: = I|Y|×|Y|.
Since W∗ satisfies L([z]TV ,:, [Y]TV ,:) = 0, we know W∗ is the optimal weight matrix that we can learn under TV . Now consider an arbitrary target node v ∈ DV with class label yv ∈ Y , and a unit structural perturbation leveraging gambit node u ∈ V with degree da that affects the predictions zv of node v made by GNN f (2)s . Without loss of generality, we assume node v has yv = 1. Note that we will only discuss the case of direct structural perturbation to the 1-hop neighborhood N(v) of target node
v, as indirect perturbations do not affect the predictions zv for node v produced by a single GNN layer. Denote ∆Ā = Ā′ − Ā as the change in the adjacency matrix Ā before and after the attack. Similar to Proof 1, for simplicity of derivation, we assume that the perturbation does not change the row-stochastic normalization of Ā, and we use δ1 to denote addition (δ1 = 1) or removal (δ1 = −1) of a homophilous edge to node v, and use δ2 to denote addition or removal of a heterophilous edge to node v.
In the case of δ1 = ±1 and δ2 = 0, we have[( (1− α)∆Ā + αI ) X ] v,: = (1− α)δ1 da|Y| [ ((|Y| − 1)p+ 1) (1− p) · · · (1− p) ]
and the change in the CM-type attack loss LCMatk before and after the perturbation can be derived as
∆LCMatk = ((1− α)|Y|+ α− 1)δ1 da(α(h− 1)− h)|Y|+ da . (19)
In the case of δ1 = 0 and δ2 = ±1, we have[( (1− α)∆Ā + αI ) X ] v,: = (1− α)δ2 da|Y| [ (1− p) (|Y| − 1)p+ 1 · · · (|Y| − 1)p+ 1 ]
and the change in the CM-type attack loss LCMatk before and after the perturbation can be derived as
∆LCMatk = (α− 1)(|Y| − 1)δ2
da(α(h− 1)− h)|Y|+ da . (20)
From Eq. (19) and Eq. (20), the change in the CM-type attack loss ∆LCMatk for GNN layer f(A,X;α) considering both δ1 and δ2 can be written as
∆LCM,fatk = ((1− α)|Y|+ α− 1)δ1 da(α(h− 1)− h)|Y|+ da + (α− 1)(|Y| − 1)δ2 da(α(h− 1)− h)|Y|+ da . (21)
Layer fs(A,X) = ĀsXW. This formulation is a special case of the previously discussed f(A,X;α) formulation when α = 11+da .
In the case of δ1 = ±1 and δ2 = 0, from Eq. (19), we have the change in the CM-type attack loss LCMatk before and after the perturbation as
∆LCMatk = − (|Y| − 1)δ1
da(h|Y| − 1) + |Y| − 1 . (22)
In the case of δ1 = 0 and δ2 = ±1, from Eq. (20), we have the change in the CM-type attack loss LCMatk before and after the perturbation as
∆LCMatk = (|Y| − 1)δ2
da(h|Y| − 1) + |Y| − 1 . (23)
From Eq. (22) and Eq. (23), the change in the CM-type attack loss ∆LCMatk for GNN layer fs(A,X) considering both δ1 and δ2 can be written as
∆LCM,fsatk = − (|Y| − 1)δ1
da(h|Y| − 1) + |Y| − 1 + (|Y| − 1)δ2 da(h|Y| − 1) + |Y| − 1
(24)
Comparison of increase in attack loss ∆LCMatk . Solving the following system of inequalities for variable α ∆LCM,fsatk > ∆L CM,f atk > 0 α, h ∈ [0, 1] |Y| ≥ 2 da, |Y| ∈ Z+ δ1, δ2 ∈ {−1, 0, 1}
(25)
we get the valid range of α as 1 da+1 < α < 1, when 0 ≤ h < 1|Y| and 0 < da < 1−|Y| h|Y|−1 and δ1 < δ2 0 ≤ α < 1da+1 , when 0 ≤ h < 1 |Y| and da > 1−|Y| h|Y|−1 and δ1 > δ2
1 da+1 < α < 1, when 1|Y| ≤ h ≤ 1 and δ1 < δ2 . (26)
From the solution in Eq. (26), we observe that when h > 1|Y| , a unit perturbation increasing Latk as discussed in Theorem 1 (i.e. δ1 = −1 and δ2 = 0, or δ1 = 0 and δ2 = 1) will satisfy the condition δ1 < δ2, and thus lead to a strictly smaller increase ∆LCM,fatk in the attack loss for layer f(A,X;α) than the increase ∆LCM,fsatk for layer fs(A,X) if α > 1 da+1 .

heading: D DETAILED EXPERIMENTAL SETUPS AND HYPERPARAMETERS
text: D.1 INSTANTIATIONS OF DESIGN ON GNNS
We explicitly demonstrate how the heterophilous design outlined in Eq. (1) of §3.2 are instantiated in various GNN models used in our experiments. In particular, we highlight how these GNN architectures allow separate aggregations of the ego- and neighbor-embeddings.
• In H2GCN, a final hidden representation is computed for each node v ∈ V through r(final)v = CONCAT(r(0)v , r (1) v , ..., r (K) v ), where r (0) v is the non-linear ego-embedding of node features and
r (k) v are the intermediate representations aggregated in the k-th layer, where k ∈ (1, ...,K). By interpreting the update rule’s CONCAT as the ENC operation, AGGR1 as the skip connection to the ego-embedding of node features, and the concatenation of the intermediate representations as AGGR2, the ego- and neighbor-embeddings are separately aggregated and the identified heterophilous design is recovered.
• Similarly, GraphSAGE (with mean aggregator) utilizes a concatenation-based encoding scheme through their update of
r(k)v = σ ( CONCAT ( r(k−1)v , MEAN ( {r(k−1)u ,∀u ∈ N(i)} )) ·W ) where AGGR1(·) = r(k−1)u , AGGR2 is the mean function and ENC(x1,x2) = σ(CONCAT(x1,x2)· W).
• GPR-GNN embeds each node feature vector separately with a fully connected layer to compute r (0) v (or H (0) v: as in the original paper), similar to H2GCN, and then updates each node’s hidden
representations through a weighted sum of all k-th hop layers around the ego-node, where k ∈ (0, 1, ...,K). By interpreting the summation as the ENC operation as well as AGGR1(·) = γ0H(0) and AGGR2(·) = ∑K k=1 γkÃ k symH
(k−1), where γ is a vector denoting the weights associated with each k-hop ego network, the aggregation of the ego- and neighbor-embeddings can be decoupled and thus GPR-GNN also satisfies the heterophilous design.
• FAGCN follows a similar update function to GPR-GNN with
h (l) i = εh (0) i + ∑ j∈N(i) αGij√ didj h (l−1) j
where h(0)i , equivalent to r (0) i in Eq. (1), represents the non-linear ego-embedding and α G ij is a proportionality constant measuring the ratio of low and high frequency components. The relationship between FAGCN and the proposed heterophilous design can similarly be inferred by interpreting the sum as the ENC operation, AGGR1(·) = εh(0)i as a weighted skip connection to the ego-embedding of feature vectors, and the weighted sum of embeddings within the neighborhood N(i) of node i ∈ V as AGGR2(·).
• CPGNN formulates the update function of belief vectors B̄(k) after the k-th propagation layer as B̄(k) = B̄(0) + AB̄(k−1)H̄, where B̄(0) consists of prior belief vectors for each node (which can be seen as the ego-embedding r(0)i in Eq. (1)), and H̄ is the linearized compatibility matrix. The heterophilous design can be recovered by letting AGGR1(·) = B̄(0) as a skip connection to the prior belief, AGGR2(·) = AB̄(k−1)H̄, and the ENC operation as the summation.

heading: D.2 MORE DETAILS ON THE EXPERIMENTAL SETUP
text: Benchmark Implementations. Our empirical framework is built on DeepRobust (Li et al., 2020b), Python Fire1 and signac2. We incorporated the following implementations of GNN models in our framework. For GNNGuard and GCN-SVD, there are some implementation ambiguities, which we discuss in the next paragraph.
H2GCN (Zhu et al., 2020) https://github.com/GemsLab/H2GCN GraphSAGE (Hamilton et al., 2017) Implemented on top of https://github.com/GemsLab/H2GCN CPGNN (Zhu et al., 2021) https://github.com/GemsLab/CPGNN GNNGuard (Zhang & Zitnik, 2020) https://github.com/mims-harvard/GNNGuard ProGNN (Jin et al., 2020b) https://github.com/ChandlerBang/Pro-GNN GCN-SVD (Entezari et al., 2020) https://github.com/DSE-MSU/DeepRobust/blob/master/examples/graph/test_gcn_svd.py GAT (Veličković et al., 2018) https://github.com/DSE-MSU/DeepRobust/blob/master/deeprobust/graph/defense/gat.py GCN (Kipf & Welling, 2017) https://github.com/DSE-MSU/DeepRobust/blob/master/deeprobust/graph/defense/gcn.py
Notes on the GNNGuard and GCN-SVD Implementations. We note that there are ambiguities in the implementations of GNNGuard (Zhang & Zitnik, 2020) and GCN-SVD (Entezari et al., 2020), which can lead to different variants with different performance and robustness, as we show in Table 6.
For GNNGuard, the ambiguity comes from different interpretations of Eq. (4) in the original paper (Zhang & Zitnik, 2020): we consider the authors’ original implementation as variant (I), and the model described in the original paper as variant (II), which we implement by building on the authors’ implementation.
1https://github.com/google/python-fire 2https://signac.io
Table 6 shows that the differences in accuracy between the two variants are in most cases less than 2%, while in many cases variant (II) shows better accuracy compared to variant (I), especially in experiments against Metattack. Thus, we use variant (II) as the default implementation for the empirical evaluations in §5.
For GCN-SVD, the ambiguity comes from the order of applying the preprocessing and low-rank approximation for the adjacency matrix A, which is not discussed in the original paper (Entezari et al., 2020).
• Variant (I): Since the original authors’ implementation is not publicly available, we consider the implementation provided in DeepRobust (Li et al., 2020b) as variant (I): it first calculates the rank-k approximation Ã of A, and then generates the preprocessed adjacency matrix Âs = D̃ −1/2 s (Ã + I)D̃ −1/2 s = D̃ −1/2 s ÃsD̃ −1/2 s , which is then processed by a GCN (Kipf & Welling,
2017). However, as the identity matrix I is added into Ã after the low-rank approximation, the diagonal elements of the resulting Âs matrix (i.e., the weights for the self-loop edges in the graph) can become significantly larger than the off-diagonal elements, especially when the rank k is low. As a result, this order of applying the preprocessing and low-rank approximation inadvertently adopts Design 1 which we identified; we have shown in Theorem 3 that even merely increasing the weights α for the ego-embedding in the linear combination ENC in Eq. (1) can lead to reduced attack loss Latk under structural perturbations. • Variant (II): In variant (II), we consider the opposite order where we first add the identity matrix I (self-loops) into the original adjacency matrix A, then we perform the low-rank approximation, and finally we symmetrically normalize the low-rank matrix Ãs to generate the preprocessed Âs used by a GCN model. This order allows the diagonal elements to be more on par in magnitude with the off-diagonal elements. As an example, on Citeseer, when using variant (I) with rank k = 5, the average magnitude of the diagonal elements of the resulting Âs can be 22.3 times the average magnitude of the off-diagonal elements; when using variant (II) instead, the average magnitude of the diagonal elements is only 9.0 times that of the off-diagonal elements.
In Table 6, we report the performance of the two variants of GCN-SVD under the experimental settings considered in §5 with rank k ∈ {5, 10, 15, 50}: Variant (I), with our first design implicitly builtin, has in most cases significantly higher performance than variant (II), especially on homophilous datasets and when rank k is low. These results further demonstrate the effectiveness of Design 1 that we identified. To enable a clear perspective of the performance and robustness improvement brought by Design 1, in our empirical analysis in §5, on top of the low-rank approximation vaccination we adopt variant (II) as the default implementation.
Attack Implementations. We incorporate the following implementations of attacks from DeepRobust (Li et al., 2020b) to our empirical framework.
NETTACK (Zügner et al., 2018) https://github.com/DSE-MSU/DeepRobust/blob/master/deeprobust/graph/targeted_attack/nettack.py Metattack (Zügner & Günnemann, 2019a) https://github.com/DSE-MSU/DeepRobust/blob/master/deeprobust/graph/global_attack/mettack.py
More Details on the Attack Setup. For NETTACK, we randomly select 60 nodes from the graph as the target nodes for each set of perturbations, instead of the GCN-based target selection approach as in (Zügner et al., 2018): the approach in (Zügner et al., 2018) only selects nodes that are correctly classified by GCN (Kipf & Welling, 2017) on clean data, thus introducing unfair advantages towards GCN, especially on heterophilous datasets where GCN can exhibit significantly inferior accuracy to models like GraphSAGE (Zhu et al., 2020). For the experiments in §5.1, we use a budget of 1 perturbation per target node to match the setups of our theorems; for the benchmark study in §5.2, we use an attack budget equal to a node’s degree and allow direct attacks on target nodes. For Metattack, we budget the attack as 20% of the number of edges in each dataset, and we use the Meta-Self variant as it shows the most destructiveness (Zügner & Günnemann, 2019a).
More Details on Randomized Smoothing Setup. Following Bojchevski et al. (2020), we similarly set the significance level α = 0.01 (i.e., the certificates hold with probability 1 − α = 0.99), using 103 samples to estimate the predictions of the smoothed classifier f(φ(s)) for input s, and another 106 samples to obtain multi-class certificates. For the randomization scheme φ, we only consider structural perturbations where with probability p+ an new edge is added, and with probability p− an existing edge is removed. We consider multiple sets of (p+, p−) in our experiments for a finer-
grained evaluation: (1) p+ = 0.001, p− = 0.4, where both addition and deletion are allowed; (2) p+ = 0.001, p− = 0, where only addition is allowed; and (3) p+ = 0, p− = 0.4, where only deletion is allowed.
Hardware Specifications. We use a workstation with a 12-core AMD Ryzen 9 3900X CPU, 64GB RAM, and a Quadro P6000 GPU with 24 GB GPU Memory.

heading: D.3 COMBINING HETEROPHILOUS DESIGN WITH LOW-RANK APPROXIMATION
text: In this section, we provide more details on how we incorporate the low-rank approximation vaccination into the formulations of H2GCN (Zhu et al., 2020) and GraphSAGE (Hamilton et al., 2017) in order to form the hybrid methods, H2GCN-SVD and GraphSAGE-SVD.
H2GCN-SVD. From (Zhu et al., 2020), each layer in the neighborhood aggregation stage of H2GCN can be algebraically formulated as
R(k) = CONCAT ( Â2R (k−1), ÂR(k−1),R(k−1) ) , (27)
where Â = D−1/2AD−1/2 is the symmetrically normalized adjacency matrix without selfloops; Â2 = D −1/2 2 A2D −1/2 2 is the symmetrically normalized 2-hop graph adjacency matrix A2 ∈ {0, 1}|V|×|V|, with [A2]u,v = 1 if v is in the 2-hop neighborhood N2(u) of node u; R(k) are the node representations after the k-th layer, and CONCAT is the column-wise concatenation function. For H2GCN-SVD, we replace Â and Â2 in Eq. (27) respectively with the low-rank approximations of Ã and Ã2, which are both postprocessed to be symmetrically normalized.
GraphSAGE-SVD. From (Hamilton et al., 2017), each layer in GraphSAGE can be algebraically formulated as R(k) = σ ( CONCAT ( ĀR(k−1),R(k−1) ) ·W(k−1) ) , (28)
where Ā is the row-stochastic graph adjacency matrix without self-loops; R(k) are the node representations after the k-th layer; CONCAT is the column-wise concatenation function; W(k) is the learnable weight matrix for the k-th layer, and σ is the non-linear activation function (ReLU). For GraphSAGE-SVD, we replace Ā in Eq. (28) with the low-rank approximation of the adjacency matrix Ã, postprocessed by row-stochastic normalization. Note that we do not enable the neighborhood sampling function for the GraphSAGE and GraphSAGE-SVD models tested in this work, as noted in Appendix D.4.
D.4 HYPERPARAMETERS
• H2GCN-SVD
Initialization Parameters:
– adj svd rank: best k chosen from {5, 50} for each dataset Training Parameters: – early stopping: Yes – train iters: 200 – patience: 100
• GraphSAGE-SVD
Initialization Parameters:
– adj nhood: [’1’] – network setup: I-T1-G-V-C1-M64-R-T2-G-V-C2-MO-R
– adj norm type: rw – adj svd rank:
best k chosen from {5, 50} for each dataset
Training Parameters: – early stopping: Yes – train iters: 200 – patience: 100
• H2GCN
Initialization Parameters: (default parameters)
Training Parameters: – early stopping: Yes – train iters: 200 – patience: 100 – lr: 0.01
• GraphSAGE
Initialization Parameters:
– adj nhood: [’1’] – network setup: I-T1-G-V-C1-M64-R-T2-G-V-C2-MO-R
– adj norm type: rw
Training Parameters: – early stopping: Yes – train iters: 200 – patience: 100 – lr: 0.01
• CPGNN
Initialization Parameters:
– network setup: M64-R-MO-E-BP2
Training Parameters: – early stopping: Yes – train iters: 400 – patience: 100 – lr: 0.01
• GPR-GNN
Initialization Parameters:
– nhid: 64 – alpha: 0.9, which is chosen from the best α ∈ {0.1, 0.2, 0.5, 0.9} on all datasets
Training Parameters: – train iters: 200 – lr: 0.01
• FAGCN
Initialization Parameters:
– nhid: 64 – alpha: 0.9, which is chosen from the best α ∈ {0.1, 0.2, 0.5, 0.9} on all datasets
– dropout: 0.5
Training Parameters: – early stopping: Yes – lr: 0.01
• GNNGuard
Initialization Parameters:
– nhid: 64 – dropout: 0.5 – base model: GCN for variant (I); GCN-fixed for variant (II) (default).
Training Parameters: – train iters: 81 – lr: 0.01
• ProGNN
Initialization Parameters: – nhid: 64 – dropout: 0.5
Training Parameters: – epochs: 400 – lr: 0.01 – lr adj: 0.01 – weight decay: 5e-4 – alpha: 5e-4 – beta: 1.5 – gamma: 1 – lambda : 0 – phi: 0 – outer steps: 1 – innter steps: 2
• GCN-SVD
Initialization Parameters:
– nhid: 64 – k: best k chosen from {5, 10, 15, 50} for
each dataset – dropout: 0.5 – svd solver: eye-svd (for variant (II) only)
Training Parameters: – train iters: 200 – weight decay: 5e-4 – lr: 0.01
• GCN
Initialization Parameters (in class MultiLayerGCN):
– nhid: 64 – nlayer: 2
Training Parameters: – train iters: 200 – lr: 0.01 – weight decay: 5e-4
• GAT
Initialization Parameters – nhid: 8 – heads: 8 – dropout: 0.5
Training Parameters: – early stopping: Yes – train iters: 1000 – patience: 100 – lr: 0.01 – weight decay: 5e-4
• MLP
Initialization Parameters: (in class H2GCN): – network setup: M64-R-D0.5-MO
Training Parameters:
– early stopping: Yes – train iters: 200 – patience: 100 – lr: 0.01

heading: D.5 DATASETS

heading: Dataset and Unidentifiability.
text: • Heterophilous Datasets: FB100 (Traud et al., 2012) is a set of 100 university friendship network snapshots from Facebook in 2005 (Lim et al., 2021), from which we use one network. Each node is labeled with the reported gender, and the features encode education and accommodation. Data is sent to the original authors (Lim et al., 2021) in an anonymized form. Though the dataset contains limited demographic (categorical) information volunteered by users on their individual Facebook pages, we manually inspect the dataset and confirm that the anonymized dataset is not recoverable and thus not identifiable. Also, no offensive content is found within the data. Snap Patents (Leskovec et al., 2005; Leskovec & Krevl, 2014) is a utility patent citation network. Node labels reflect the time the patent was granted, and the features are derived from the patent’s metadata. The dataset is maintained by the National Bureau of Economic Research, and is freely available for download3. Neither personally identifiable information nor offensive content is identified when we manually inspect the dataset.
• Homophilous Datasets: Cora (McCallum, 2000) and Citeseer (Sen et al., 2008) datasets are scientific publication citation networks, whose labels categorize the research field, and features indicate the absence or presence of the corresponding word from the dictionary. No personally identifiable information or offensive content is identified when we manually inspect both datasets.
Downsampling. For better computational tractability, we sample a subset of the Snap Patents data using a snowball sampling approach (Goodman, 1961), where a random 20% of the neighbors for each traversed node are kept. We provide the pseudocode for the downsampling process in Algorithm 1.
Algorithm 1: Downsampling Algorithm For Snap Patents Input: Graph to sample G
Number of nodes to sample N Sampling ratio p
Output: Downsampled graph G′ 1 initialization /* Initialize a queue bfsquene for Breadth First Search, and a list nodessampled for storing sampled nodes */ 2 bfsquene ← QUEUE() 3 nodessampled ← LIST() /* Start BFS with a random node from the largest connected component
in G; RANDOM(array, n) returns n elements from an array with equal probability without replacement */
4 nodestarting ← RANDOM(LARGESTCONNECTEDCOMPONENT(G), 1) 5 push nodestarting into bfsquene 6 while LENGTH(nodessampled) ¡ N do 7 node← bfsquene.pop() 8 neighbors← one hop neighbors of node 9 neighborsdrawn ← RANDOM(neighbors, p× LENGTH(neighbors))
10 for neighbor ∈ neighborsdrawn do 11 if neighbor /∈ nodessampled then 12 append neighbor to nodessampled 13 push neighbor into bfsquene 14 end 15 end 16 end 17 G′← subgraph induced by nodessampled 18 return G′
3https://www.nber.org/research/data/us-patents

heading: E.1 DETAILED RESULTS FOR EVALUATION ON EMPIRICAL ROBUSTNESS

heading: E DETAILED EXPERIMENT RESULTS

heading: E.2 DETAILED RESULTS FOR EVALUATION ON CERTIFIABLE ROBUSTNESS

heading: E.3 COMPARISON BETWEEN CERTIFIABLE AND EMPIRICAL ROBUSTNESS
text: While the evaluations on both empirical and certifiable robustness have shown that methods featuring the design have shown largely improved robustness compared to the best-performing unvaccinated method, we find that the robustness rankings for methods under certifiable and empirical robustness are different; previous results from Geisler et al. (2020) have also shown discrepancy in certifiable and empirical robustness rankings.
We think this discrepancy may be attributed to multiple factors. Firstly, the radius on which certificates can be issued with randomized smoothing may not cover the radius of perturbations we allowed for NETTACK and Metattack in §5.2, which are more than tens of edges for Metattack and high-degree nodes in NETTACK (where we use an attack budget equal to the degree of the target node). Furthermore, even for low-degree nodes, attacks are much more inclined to introduce new edges instead of removing existing ones, as we have shown in §5.1, which can make it much harder to obtain certificates (Bojchevski et al., 2020). In our case, the methods we evaluated generally have r̄a < 1, meaning that for most nodes there are no certificates to cover even the addition of a single edge. Thus, methods which display higher certifiable robustness under smaller perturbations may not keep their robustness under larger perturbations by empirical attacks. Moreover, while we are evaluating on randomized smoothed models f(φ(s)) to measure certifiable robustness, in empirical
robustness we are evaluating on the robustness of the base models f(s), which may differ from the robustness of the randomized smoothed models. Lastly, it is worth noting that the lack of certification for a model within certain radius does not imply a vulnerability to all adversarial attacks within that radius; the model may still be robust against many attacks within that radius. Similarly, it is also possible for existing certification approaches to miss some certifiable cases. Taking all these factors into account, we believe that while evaluations on certifiable robustness provide complementary perspectives to evaluations empirical robustness, at the current stage it cannot replace the evaluations on empirical robustness, and the relation between certifiable and empirical robustness remains as a question for future works.

heading: E.4 COMPLEXITY AND RUNTIME OF HETEROPHILOUS VACCINATION
text: Another benefit of adopting heterophily-inspired design for boosting robustness of GNNs is their smaller computational overhead compared to existing vaccination mechanisms, especially vaccinations based on low-rank approximation. As our identified design can be applied as simple architectural changes on top of an existing GNN, they usually maintain the same order of computational complexity as the base model. For example, adding the heterophilous design to GCN (Kipf & Welling, 2017) results in an architecture similar to GraphSAGE (Hamilton et al., 2017); both have the same order of computational complexity as O(|V| + |E|) by leveraging the sparse connectivity of most real-world graphs. Low-rank approximation-based vaccination, on the other hand, approximates the adjacency matrix of a graph by an SVD, resulting in an adjusted low-rank adjacency matrix Ã based on which the GNN runs. However, not only is computing an SVD potentially costly (O(|V|3) in general), but in most cases it also results in a dense Ã (in contrast to the sparse original adjacency matrix), thus increasing the complexity of each iteration of the GNN.
Table 10 shows the runtime of 200 training iterations of each model. We observe that models with the heterophilous design have the smallest runtime among all vaccinated models. Even for methods based on the same implementation, H2GCN and GraphSAGE are still 3-4 times faster than the corresponding H2GCN-SVD and GraphSAGE-SVD methods. For fair runtime measurements, we measure the runtime of each model on an Amazon EC2 instance with instance type as p3.2xlarge, which features an 8-core CPU, 61 GB Memory, and a Tesla V100 GPU with 16 GB GPU Memory.
E.5 ADDITIONAL RESULTS ON STRUCTURAL ATTACKS ON HETEROPHILOUS GRAPHS

Review summaries:
This paper investigates the change of the degree of homophily under adversarial perturbation of graphs from both theoretical and empirical points of view. In theory, this paper proves under simplified settings, for homophily graphs, the attack will increase the level of heterophily, and for heterophily graphs, this change will depend on the degree. In practice, this paper advocates a design principle of GNNs that separate aggregators for ego- and neighbor-embeddings.

This paper studies the relationship between graph heterophily and robustness of GNNs. The authors claim that effective structural attacks on GNNs for homophilious graphs lead to increased heterophily level, while for heterophilious graphs they alter the homophily level contingent on node degrees. The paper shows that a heterophily-inspired design for GNN, separate aggregators for ego- and neighbor-embeddings, could improve the robustness of several GNN architectures.

The authors present a study on the robustness on GNNs on homophilic and heterophilic datasets. The authors first show in a stylized theoretical setting that GNNs with a *heterophilic design* (separate aggregation of ego- and neighborhood representations) should show increased robustness against structure perturbations. Then, they follow up with an experimental study where they compare homophily-based GNNs with GNNs following a heterophilic design on both homophilic and heterophilic datasets. In addition, the authors augment three GNN models with the SVD preprocessing by [Entezari et al. 2020] in order to improve robustness. The authors study both empirical as well as certifiable robustness (via randomized smoothing). In general, GNNs with heterophilic design show higher empirical and certifiable robustness across datasets.