Title: SUPERVISED COMMUNITY DETECTION WITH LINE GRAPH NEURAL NETWORKS

Abstract: Community detection in graphs can be solved via spectral methods or posterior inference under certain probabilistic graphical models. Focusing on random graph families such as the stochastic block model, recent research has unified both approaches and identified both statistical and computational detection thresholds in terms of the signal-to-noise ratio. By recasting community detection as a node-wise classification problem on graphs, we can also study it from a learning perspective. We present a novel family of Graph Neural Networks (GNNs) for solving community detection problems in a supervised learning setting. We show that, in a data-driven manner and without access to the underlying generative models, they can match or even surpass the performance of the belief propagation algorithm on binary and multiclass stochastic block models, which is believed to reach the computational threshold in these cases. In particular, we propose to augment GNNs with the non-backtracking operator defined on the line graph of edge adjacencies. The GNNs are achieved good performance on real-world datasets. In addition, we perform the first analysis of the optimization landscape of using (linear) GNNs to solve community detection problems, demonstrating that under certain simplifications and assumptions, the loss value at any local minimum is close to the loss value at the global minimum/minima.

1. INTRODUCTION:

Graph inference problems encompass a large class of tasks and domains, from posterior inference in probabilistic graphical models to community detection and ranking in generic networks, image segmentation or compressed sensing on non-Euclidean domains. They are motivated both by practical applications, such as in the case of PageRank, and also by fundamental questions on the algorithmic hardness of solving such tasks.

From a data-driven perspective, these problems can be formulated in unsupervised, semi-supervised or supervised learning settings. In the supervised case, one assumes a dataset of graphs with labels on their nodes, edges or the entire graphs, and attempts to perform node-wise, edge-wise and graph-wise classification by optimizing a loss over a certain parametric class, e.g. neural networks. Graph Neural Networks (GNNs) are natural extensions of Convolutional Neural Networks to graph-structured data, and have emerged as a powerful class of algorithms to perform complex graph inference leveraging labeled data (Gori et al., 2005; Bronstein et al., 2017) (and references therein). In essence, these neural networks learn cascaded linear combinations of intrinsic graph operators interleaved with node-wise (or edge-wise) activation functions. Since they utilize intrinsic graph operators, they can be applied to varying input graphs, and they offer the same parameter sharing advantages as their CNN counterparts.

In this work, we focus on community detection problems, a wide class of node classification tasks that attempt to discover a clustered, segmented structure within a graph. The algorithmic approaches to this problem include a rich class of spectral methods, which take advantage of the spectrum of certain operators defined on the graph, as well as approximate message-passing methods such as belief propagation (BP), which performs approximate posterior inference under predefined graphical models (Decelle et al., 2011) . Focusing on the supervised setting, we study the ability of GNNs to approximate, generalize or even improve upon these class of algorithms. Our motivation is two-fold. On the one hand, this problem exhibits algorithmic hardness on some settings, opening up the possibility to discover more efficient algorithms than the current ones. On the other hand, many practical scenarios fall beyond pre-specified probabilistic models, requiring data-driven solutions.

We propose modifications to the GNN architecture, which allow it to exploit edge adjacency information, by incorporating the non-backtracking operator of the graph. This operator is defined over the edges of the graph and allows a directed flow of information even when the original graph is undirected. It was introduced to community detection problems by Krzakala et al. (2013) , who propose a spectral method based on the non-backtracking operator. We refer to the resulting GNN model as a Line Graph Neural Network (LGNN). Focusing on important random graph families exhibiting community structure, such as the stochastic block model (SBM) and the geometric block model (GBM), we demonstrate improvements in the performance by our GNN and LGNN models compared to other methods, including BP, even in regimes within the so-called computational-to-statistical gap. A perhaps surprising aspect is that these gains can be obtained even with linear LGNNs, which become parametric versions of power iteration algorithms.

We want to mention that besides community detection tasks, GNN and LGNN can be applied to other node-wise classification problems too. The reason we are focusing on community detection problems is that this is a relatively well-studied setup, for which different algorithms have been proposed and where computational and statistical thresholds have been studied in several scenarios. Moreover, synthetic datasets can be easily generated for community detection tasks. Therefore, we think it is a nice setup for comparing different algorithms, besides its practical values.

The good performances of GNN and LGNN motivate our second main contribution: the analysis of the optimization landscape of simplified and linear GNN models when trained with planted solutions of a given graph distribution. Under reparametrization, we provide an upper bound on the energy gap controlling the energy difference between local and global minima (or minimum). With some assumptions on the spectral concentration of certain random matrices, this energy gap will shrink as the size of the input graphs increases, which would mean that the optimization landscape is benign on large enough graphs.

Summary of Main Contributions::

• We propose an extension of GNNs that operate on the line graph using the non-backtracking operator, which yields improvements on hard community detection regimes.

• We show that on the stochastic block model we reach detection thresholds in a purely data-driven fashion, in the sense that our results improve upon belief propagation in hard SBM detection regimes, as well as in the geometric block model.

• We perform the first analysis of the learning landscape of GNN models, showing that under certain simplifications and assumptions, they exhibit a form of "energy gap", where local mimima are confined in low-energy configurations.

• We show that our model can perform well on community detection problems with real-world datasets.

2. PROBLEM SETUP:

We are interested in a specific class of node-classification tasks in which given an input graph G = (V, E), a labeling y : V → {1, . . . , C} that encodes a partition of V into C communities is to be predicted at each node. We assume that a training set {(G t , y t )} t≤T is given, with which we train a model that predicts ŷ = Φ(G, θ) by minimizing

L(θ) = 1 T t≤T (Φ(G t , θ), y t )

Since y encodes a partition of C groups, the specific label of each node is only important up to a global permutation of {1, . . . , C}. Section 4.3 describes how to construct loss functions with such a property. A permutation of the observed nodes translates into the same permutation applied to the labels, which justifies models Φ that are equivariant to permutations. Also, we are interested in inferring properties of community detection algorithms that do not depend on the specific size of the graphs1 . We therefore require that the model Φ accepts graphs of variable size for the same set of parameters, similar to sequential RNN or spatial CNN models.

3. RELATED WORK:

GNN was first proposed in Gori et al. (2005) ; Scarselli et al. (2009) . Bruna et al. (2013) generalize convolutional neural networks on general undirected graphs by using the graph Laplacian's eigenbasis. This was the first time the Laplacian operator was used in a neural network architecture to perform classification on graph inputs. Defferrard et al. (2016) consider a symmetric Laplacian generator to define a multiscale GNN architecture, demonstrated on classification tasks. Similarly, Kipf & Welling (2016) use a similar generator as effective embedding mechanisms for graph signals and applies it to semi-supervised tasks. This is the closest application of GNNs to our current contribution. However, we highlight that semi-supervised learning requires bootstrapping the estimation with a subset of labeled nodes, and is mainly interested in generalization within a single, fixed graph. In comparison, our setup considers community detection across a distribution of input graphs and assumes no initial labeling on the graphs in the test dataset except for the adjacency information.

There have been several extensions of GNNs by modifying their non-linear activation functions, parameter sharing strategies, and choice of graph operators (Li et al., 2015; Sukhbaatar et al., 2016; Duvenaud et al., 2015; Niepert et al., 2016) . In particular, Gilmer et al. (2017) interpret the GNN architecture as learning an approximate message-passing algorithm, which extends the learning of hidden representations to graph edges in addition to graph nodes. Recently, Velickovic et al. (2017) relate adjacency learning with attention mechanisms, and Vaswani et al. (2017) propose a similar architecture in the context of machine translation. Another recent and related piece of work is by Kondor et al. (2018) , who propose a generalization of GNN that captures high-order node interactions through covariant tensor algebra. Our approach to extend the expressive power of GNN using the line graph may be seen as an alternative to capture such high-order interactions.

Our energy landscape analysis is related to the recent paper by Shamir (2018) , which establishes an energy bound on the local minima arising in the optimization of ResNets. In our case, we exploit the properties of the community detection problem to produce an energy bound that depends on the concentration of certain random matrices, which one may hope for as the size of the input graphs increases. Finally, Zhang (2016)'s work on data regularization for clustering and rank estimation is also motivated by the success of using Bethe-Hessian-like perturbations to improve spectral methods on sparse networks. It finds good perturbations via matrix perturbations and also has successes on the stochastic block model. Yang & Leskovec (2012a) curate benchmark datasets for community detection and quantify the quality of these datasets, while Yang & Leskovec (2012b) develop new algorithms for community detection by fitting to the networks the Affliation Graph Model (AGM), a generative model for graphs with overlapping communities.

4. LINE GRAPH NEURAL NETWORKS:

This section introduces our GNN architectures that include the power graph adjacency (Section 4.1) and its extension to line graphs using the non-backtracking operator (Section 4.2), as well as the design of losses invariant to global label permutations (Section 4.3).

4.1. GRAPH NEURAL NETWORKS USING A FAMILY OF GRAPH OPERATORS:

The Graph Neural Network (GNN), introduced in Scarselli et al. (2009) and later simplified in Li et al. (2015) ; Duvenaud et al. (2015) ; Sukhbaatar et al. (2016) , is a flexible neural network architecture based on local operators on a graph G = (V, E). Given a state vector x ∈ R |V |×b on the vertices of Figure 1 . Overview of the architecture of LGNN (Section 4.2). Given a graph G, we construct its line graph L(G) with the non-backtracking operator (Figure 2 ). In every layer, the states of all nodes in G and L(G) are updated according to (2). The final states of nodes in G are used to predict node-wise labels, and the trainining is performed end-to-end using standard backpropagation with a label permutation invariant loss (Section 4.3).

G, we consider intrinsic linear operators of the graph that act locally on x, which can be represented as |V |-by-|V | matrices. For example, the adjacency matrix A is defined entry-wise by

A i1i2 = 1 if (i 1 , i 2 ) ∈ E and A i1i2 = 0 if (i 1 , i 2 ) / ∈ E, for every pair (i 1 , i 2 ) ∈ V × V .

The degree matrix D is then defined as diag(A1), i.e., D is a diagonal matrix with D ii being the number of edges that the ith node has. We can also define power graph adjacency matrices as A (j) = min(1, A2 j ), which encodes 2 j -hop neighborhoods into a binary graph. Finally, there is also the identity matrix, I. Given such a family of operators for each graph, F J A = {I, D, A, A (2) , ..., A (J) }, we define a GNN layer that maps x (k) ∈ R |V |×b k to x (k+1) ∈ R |V |×b k+1 as follows. First, we compute

z (k+1) = ρ   Oi∈F J A O i x (k) θ i   , z (k+1) = Oi∈F J A O i x (k) θ i (1)

where

θ j ∈ R b k × b k+1

2 are trainable parameters and ρ(•) is a point-wise nonlinear activation function, chosen in this work to be the ReLU function, i.e. ρ(z) = max(0, z) for z ∈ R. Then we define x (k+1) = [z (k+1) , z (k+1) ] ∈ R |V |×b k+1 as the concatenation of z (k+1) and z (k+1) . The layer thus includes linear "residual connections" (He et al., 2016) via z (k) , both to ease with the optimization when using large number of layers and to increase the expressivity of the model by enabling it to perform power iterations. Since the spectral radius of the learned linear operators in (1) can grow as the optimization progresses, the cascade of GNN layers can become unstable to training. In order to mitigate this effect, we consider spatial batch normalization (Ioffe & Szegedy, 2015) at each layer. 2 In our experiments, the initial states are set to be the degrees of the nodes, i.e., x (0) = deg(A) := A1. In addition, the model Φ(G, x (0) ) = x (K) satisfies the permutation equivariance property required for community detection: Given a permutation π among the nodes in the graph, Φ(G π , Πx (0) ) = ΠΦ(G, x (0) ), where Π is the permutation matrix associated with π.

Analogy between GNN and power iterations In our setup, spatial batch normalization not only prevents gradient blowup, but also performs the orthogonalisation relative to the constant vector, which reinforces the analogy with the spectral methods for community detection, some background of which is described in B.1. In essence, in certain regimes, the eigenvector of A corresponding to its second largest eigenvalue and the eigenvector of the Laplacian matrix, L = D -A, corresponding to its second smallest eigenvalue (i.e. the Fiedler vector), are both correlated with the community structure of the graph. Thus, spectral methods for community detection performs power iterations on these matrices to obtain the eigenvectors of interest and predicts the community structure based on them. For example, to extract the Fiedler vector of a matrix M , whose eigenvector corresponding to the smallest eigenvalue is known to be v, one performing power iterations on M = M I -M by y

(n+1) = M x (n) , x (n+1) = y (n+1) -v T vy (n+1) y (n+1) -v T vy (n+1) .

If v is a constant vector, which is the case for L, then the normalization above is precisely performed within the spatial batch normalization step.

By incorporating a family of operators into the neural network framework, the GNN can not only approximate but also go beyond power iterations. As explained in Section B.1, the Krylov subspace generated by the graph Laplacian (Defferrard et al., 2016) is not sufficient to operate well in the sparse regime, as opposed to the space generated by {I, D, A}. The expressive power of each layer is further increased by adding multiscale versions of A, although this benefit comes at the cost of computational efficiency, especially in the sparse regime. The network depth is chosen to be of the order of the graph diameter, so that all nodes obtain information from the entire graph. In sparse graphs with small diameter, this architecture offers excellent scalability and computational complexity. Indeed, in many social networks diameters are constant (due to hubs), or log(|V |), as in the stochastic block model in the constant average degree regime (Riordan & Wormald, 2010) . This results in a model with computational complexity on the order of |V | log(|V |), making it amenable to large-scale graphs.

4.2. LGNN: GNN ON LINE GRAPHS WITH THE NON-BACKTRACKING OPERATOR:

For graphs with few cycles, posterior inference can be remarkably approximated by loopy belief propagation (Yedidia et al., 2003) . As described in Section B.2, the message-passing rules are defined over the edge adjacency graph (see equation 57). Although its second-order approximation around the critical point can be efficiently approximated with a power method over the original graph, a datadriven version of BP requires accounting for the non-backtracking structure of the message-passing. In this section we describe an upgraded GNN model that exploits the non-backtracking structure. 

B (i→j),(i →j ) = 1 if j = i and j = i , 0 otherwise.

This operator enables the directed propagation of information through on the line graph and was first proposed in the context of community detection on sparse graphs in Krzakala et al. (2013) . The message-passing rules of BP can be expressed as a diffusion in the line graph L(G) using this non-backtracking operator, with specific choices of activation function that turn product of beliefs into sums.

A natural extension of the GNN architecture presented in Section 4.1 is thus to consider a second GNN defined on L(G), where B and D B = diag(B1) play the role of the adjacency and the degree matrices, respectively. This effectively defines edge features that are updated according to the edge adjacency of G. Edge and node features communicate at each layer using the edge indicator matrices P m , P d ∈ {0, 1} |V |×2|E| , defined as P mi,(i→j) = 1, P dj,(i→j) = 1, P di,(i→j) = 1, P dj,(i→j) = -1 and 0 otherwise. With the skip linear connections defined similarly, the resulting model becomes

EQUATION

where

F A = {I, D, A, A (2) , . . . , A (J) }, F B = {I B , D B , B, B (2) , . . . , B (J) }, F AB = {P m , P d },

and the trainable parameters are θ i , θ i , θ i ∈ R b k ×b k+1 and θ i ∈ R b k+1 ×b k+1 . We call such a model a Line Graph Neural Network (LGNN).

In our experiments, we set x (0) = deg(A) and y (0) = deg(B). For graph families with constant average degree d (as |V | grows), the line graph has size 2|E| ∼ O(d|V |), and is therefore feasible from the computational point of view. Furthermore, the construction of line graphs can be iterated to generate L(L(G)), L(L(L(G))), etc. to yield a graph hierarchy, which could capture high-order interactions among nodes of G. Such an hierarchical construction is related to other recent efforts to generalize GNNs (Kondor et al., 2018) .

Relationship between:

LGNN and edge feature learning approaches Several authors have proposed combining node and edge feature learning. Battaglia et al. (2016) introduce edge features over directed and typed graphs, but does not discuss the undirected case. Kearnes et al. (2016) ; Gilmer et al. (2017) learn edge features on undirected graphs using f e = g(x(i), x(j)) for an edge e = (i, j), where g is commutative on its arguments. Finally, Velickovic et al. ( 2017) learns directed edge features on undirected graphs using stochastic matrices as adjacencies (which are either row or column-normalized). However, we are not aware of works that consider the edge adjacency structure provided by the non-backtracking matrix on the line graph. With non-backtracking matrix, our LGNN can be interpreted as learning directed edge features from an undirected graph. Indeed, if each node i contains two distinct sets of features x s (i) and x r (i), the non-backtracking operator constructs edge features from node features while preserving orientation: For an edge e = (i, j), our model is equivalent to constructing oriented edge features f i→j = g(x s (i), x r (j)) and f j→i = g(x r (i), x s (j)) (where g is trainable and not necessarily commutative on its arguments) that are subsequently propagated through the graph. Constructing such local oriented structure is shown to be important for improving performance in the next section.

For comparison, we also define a linear LGNN (LGNN-L) as the the LGNN that drops the nonlinear activation functions ρ in (2), and a symmetric LGNN (LGNN-S) as the LGNN whose line graph is defined on the undirected edges of the original graph: In LGNN-S, two edges of G are connected in the line graph if and only if they share one common node; also, F = {P }, with P ∈ R |V |×|E| defined as P i,(j→k) = 1 if i = j or k and 0 otherwise.

4.3. A LOSS FUNCTION INVARIANT UNER LABEL PERMUTATION:

Let C = {1, . . . , C} denote the set of all community labels, and consider first the case where communities do not overlap. By applying the softmax function at the end, we interpret the cth dimension of the output of the models at node i as the conditional probability that the node belongs to community c, o i,c = p(y i = c |θ, G). Let G = (V, E) be the input graph and let y i be the ground truth community label of node i. Since the community structure is defined up to global permutations of the labels, we can define a loss function with respect to a given graph instance as

EQUATION

where S C denotes the permutation group of C elements. This is essentially taking the the cross entropy loss minimized over all possible permutations of C. In our experiments, we considered examples with small numbers of communities such as 2 and 5. In general scenarios where C is much larger, the evaluation of the loss function (3) can be impractical due to the minimization over S C . A possible solution is to randomly partition C labels into C groups, and then to marginalize the model outputs Finally, if communities may overlap, we can enlarge C to include subsets of communities and define the permutation group accordingly. For example, if there are two overlapping communities, we let C = {{1}, {2}, {1, 2}}, and only allow the permutation between 1 and 2 when computing the loss function as well as the overlap to be introduced in Section 6.

o i,c , c ≤ C into ōi,c = c∈c o i,c , c ∈ C,

5. ENERGY LANDSCAPE OF LINEAR GNN OPTIMIZATION:

As described in the numerical experiments, we found that the GNN models without nonlinear activations already provide substantial gains relative to baseline (non-trainable) algorithms. This section studies the optimization landscape of linear GNNs. Despite defining a non-convex objective, we prove that the landscape is "benign" under certain further simplifications, in the sense that the local minima are confined in sublevel sets of low energy.

For simplicity, we consider only the binary c = 2 case where we replace the node-wise binary cross-entropy loss by the squared cosine distance3 , assume a single feature map (b k = 1 for all k), and focus on the GNN described in Section 4.1 (although our analysis carries equally to describe the line graph version; see remarks below). We also make the simplifying assumption to replace the layer-wise spatial batch normalization by a simpler projection onto the unit 2 ball (thus we do not remove the mean). Without loss of generality, assume that the input graph G has size n, and denote by F = {A 1 , . . . , A Q } the family of graph operators appearing in (1). Each layer thus applies an arbitrary polynomial

Q q=1 θ (k)

q A q to the incoming node feature vector x (k) . Given an input node vector w ∈ R n , the network output can thus be written as

EQUATION

We highlight that this linear GNN setup is fundamentally different from the linear fully-connected neural networks (that is, neural networks with linear activation function), whose landscape has been analyzed in Kawaguchi (2016) . First, the output of the GNN is on the unit sphere, which has a different geometry. Next, since the operators in F depend on the input graph, they introduce fluctuations in the landscape. In general, the operators in F are not commutative, but by considering the generalized Krylov subspace generated by powers of F,

F K = {O 1 = A K 1 , O 2 = A 1 A K-1 2 , O 3 = A 1 A 2 A K-2 1 , . . . O Q K = A K Q }, one can reparametrize (4) as e = Q K j=1 β j O j w with β ∈ R M , with M = Q K .

Given the target y ∈ R n , the loss incurred by each pair (G, y) becomes 1 -| e,y | 2 e 2 , and therefore the population loss, when expressed in terms of β, equals

EQUATION

Y n = z n z n ∈ R M ×M , (z n ) j = O j w, y and X n = U n U n ∈ R M ×M , U n =   (O 1 w) . . . (O M w)   .

The landscape is thus specified by a pair of random matrices Y n , X n ∈ R M ×M .

Assuming that EX n 0, we write the Cholesky decomposition of EX n as EX n = R n R T n , and define

A n = R -1 n Y n (R -1 n ) T , Ān = EA n = R -1 n EY n (R -1 n ) T , B n = R -1 n X n (R -1 n ) T

, and ∆B n = B n -I n . Given a symmetric matrix K ∈ R M ×M , we let λ 1 (K), λ 2 (K), ..., λ M (K) denote the eigenvalues of K in nondecreasing order. Then, the following theorem establishes that under appropriate assumptions, the concentration of relevant random matrices around their mean controls the energy gaps between local and global minima of L. Theorem 5.1. For a given n, let

η n = (λ 1 ( Ān ) -λ 2 ( Ān )) -1 , µ n = E[|λ 1 (A n )| 6 ], ν n = E[|λ 1 (B n )| -6 ], δ n = E[ ∆B n 6 ],

and assume that all four quantities are finite. Then if β l ∈ S M -1 is a local minimum of L n , and β g ∈ S M -1 is a global minimum of L n , we have L n (β l ) ≤ (1 + ηn,µn,νn,δn ) • L n (β g ), where ηn,µn,νn,δn = O(δ n ) for given η n , µ n , ν n as δ n → 0 and its formula is given in the appendix. Corollary 5.2. If (η n ) n∈N * , (µ n ) n∈N * , (ν n ) n∈N * are all bounded sequences, and

lim n→∞ δ n = 0, then ∀ > 0, ∃n such that ∀n > n , |L n (β l ) -L n (β g )| ≤ • L n (β g ).

The main strategy of the proof is to consider the actual loss function L n as a perturbation of

Ln (β) = 1 -E Xn,Yn β T Ynβ β T EXnβ = 1 -β T EYnβ

β T EXnβ , which has a landscape that is easier to analyze and does not have poor local minima, since it is equivalent to a quadratic form defined over the sphere S M -1 . Applying this theorem requires estimating spectral fluctuations of the pair X n , Y n , which in turn involve the spectrum of the C * algebras generated by the non-commutative family F. For example, for stochastic block models, it is an open problem how the bound behaves as a function of the parameters p and q. Another interesting question is to understand how the asymptotics of our landscape analysis relate to the hardness of estimation as a function of the signal-to-noise ratio.

Finally, another open question is to what extent our result could be extended to the non-linear residual GNN case, perhaps leveraging ideas from Shamir (2018) .

6. EXPERIMENTS:

We present experiments on community detection in synthetic datasets (Sections 6.1, 6.2 and Appendix C.1) as well as real-world datasets (Section 6.3). In the synthetic experiments, the performance is measured by the overlap between predicted (ŷ) and true labels (y), which quantifies how much better than random guessing a predicted labeling is, given by 1

n u δ y(u),ŷ(u) -1 C /(1 -1 C )

, where δ is the Kronecker delta function, and this quantity is maximized over global permutations within a graph of the set of labels. In the real-world datasets, as the communies are overlapping and not balanced, the prediction accuracy is measured by 1 n u δ y(u),ŷ(u) , and the set of permutations to be maximized over is described in Section 4.3. We used Adamax (Kingma & Ba, 2014) with learning rate 0.004 across all experiments. All the neural network models have 30 layers and 8 features in the middle layers (i.e., b k = 8) for experiments in Sections 6.1 and 6.2, and 20 layers and 6 features for Section 6.3. GNNs and LGNNs have J = 2 across the experiments except the ablation experiments in Section C.3.4 

6.1. BINARY STOCHASTIC BLOCK MODEL:

The stochastic block model is a random graph model with planted community structure. In its simplest form, the graph consists of |V | = n nodes, which are partitioned into C communities, that is, each node is assigned a label y ∈ {1, ..., C}. An edge connecting any two vertices u, v is drawn independently at random with probability p if y(v) = y(u), and with probability q otherwise. In the binary case (i.e. C = 2), the sparse regime, where p, q 1/n, is well understood and provides an initial platform to compare the GNN and LGNN with provably optimal recovery algorithms (Appendix B). We consider two learning scenarios. In the first scenario, we choose different pairs of p and q, and train the models for each pair separately. In particular, for each pair of (p i , q i ), we sample 6000 graphs under G ∼ SBM (n = 1000, p i , q i , C = 2) and then train the models for each i. In the second scenario, reported in Appendix C.2, we train a single set of parameters θ from a set of 6000 graphs sampled from a mixture of SBM with different pairs of (p i , q i ), and average degree. Importantly, his setup shows that our models are not simply approximating known algorithms such as BP for particular SBM parameters, since the parameters vary in this dataset.

For the first scenario, we chose five different pairs of (p i , q i ) while fixing p i + q i , thereby corresponding to different signal-to-noise ratios (SNRs). Figure 3 reports the performance of our models on the binary SBM model for the different SNRs, compared with baseline methods including BP, spectral methods using the normalized Laplacian and the Bethe Hessian as well as Graph Attention Networks (GAT)5 from Velickovic et al. (2017) . We observe that both GNN and LGNN reach the performance of BP. In addition, even the linear LGNN achieves a performance that is quite close to that of BP, in accordance to the spectral approximations of BP given by the Bethe Hessian (see supplementary), and significantly outperforms performing 30 power iterations on the Bethe Hessian or the normalized Laplacian, as was done in the spectral methods. We also notice that our models outperform GAT in this task. We ran experiments in the dissociative case (q > p), as well as with C = 3 communities and obtained similar results, which are not reported here. Table 1 : Performance of different models on 5-community dissociative SBM graphs with n = 400, C = 5, p = 0, q = 18/n, corresponding to an average degree of 14.5. The first row gives the average overlap across test graphs, and the second row gives the graph-wise standard deviation of the overlap.

6.2. COMPUTATIONAL-TO-STATISTICAL THRESHOLDS IN THE SBM:

In SBM with fewer than 4 communities, it is known that BP provably reaches the information-theoretic threshold (Abbe, 2017; Massoulié, 2014; Coja-Oghlan et al., 2016) . The situation is different for k > 4, where it is conjectured that a gap emerges between the theoretical performance of MLE estimators and the performance of any polynomial-time estimation procedure (Decelle et al., 2011) .

In this context, one can use the GNN models to search the space of the generalizations of BP, and attempt to improve upon the detection performance of BP for scenarios where the SNR falls within the computational-to-statistical gap. Table 1 presents results for the 5-community dissociative SBM, with n = 400, p = 0 and q = 18/n. The SNR in this setup is above the information-theoretic threshold but below the asymptotic threshold above which BP is able to detect (Decelle et al., 2011) . Note that since p = 0, this also amounts to a graph coloring problem.

We see that the GNN and LGNN models outperform BP in this experiment, indeed opening up the possibility to reduce the computation-information gap. That said, our model may taking advantage of finite-size effects, which will vanish as n → ∞. The asymptotic study of these gains is left for future work. In terms of average test accuracy, LGNN has the best performance. In particular, it outperforms the symmetric version of LGNN, emphasizing the importance of the non-backtracking matrix used in LGNN. Although equipped with the attention mechanism, GAT does not explicitly incorporate in itself the degree matrix, the power graph adjacency matrices or the line graph structure, and has inferior performance compared with the GNN and LGNN models. Further ablation studies on GNN and LGNN are described in Section C.3.

6.3. REAL DATASETS FROM SNAP:

We now compare the models on the SNAP datasets, whose domains range from social networks to hierarchical co-purchasing networks. We obtain the training set as follows. For each SNAP dataset, we start by focusing only on the 5000 top quality communities provided by the dataset. We then identify edges (i, j) that cross at least two different communities. For each of such edges, we consider pairs of communities C 1 , C 2 such that i / ∈ C 2 and j / ∈ C 1 , i ∈ C 1 , j ∈ C 2 , and extract the subset of nodes determined by C 1 ∪ C 2 together with the edges among them. The resulting graph is connected since each community is connected. Finally, we divide the dataset into training and testing sets by enforcing that no community belongs to both the training and the testing set. In our experiment, due to computational limitations, we restrict our attention to the three smallest datasets in the SNAP collection (Youtube, DBLP and Amazon), and we restrict the largest community size to 200 nodes, which is a conservative bound.

We compare the performance of GNN and LGNN models with GAT as well as the Community-Affiliation Graph Model (AGM), which is a generative model proposed in Yang & Leskovec (2012b) that captures the overlapping structure of real-world networks. Community detection can be achieved by fitting AGM to a given network, which was shown to outperform some state-of-the-art algorithms. Table 2 compares the performance, measured with a 3-class (C = {{1}, {2}, {1, 2}}) classification accuracy up to global permutation 1 ↔ 2. GNN, LGNN, LGNN-S and GAT yield similar results and outperform AGMfit. It further illustrates the benefits of data-driven models that strike the right balance between expressivity and structural design. 

7. CONCLUSION:

In this work, we have studied data-driven approaches to supervised community detection with graph neural networks. Our models achieve comparable performance to BP in binary SBM for various SNRs, and outperform BP in the sparse regime of 5-class SBM that falls between the computationalto-statistical gap. This is made possible by considering a family of graph operators including the power graph adjacency matrices, and importantly by introducing the line graph equipped with the non-backtracking matrix. We also provided a theoretical analysis of the optimization landscapes of simplified linear GNN for community detection and showed the gap between the loss value at local and global minima are bounded by quantities related to the concentration of certain random matricies.

One word of caution is that our empirical results are inherently non-asymptotic. Whereas models trained for given graph sizes can be used for inference on arbitrarily sized graphs (owing to the parameter sharing of GNNs), further work is needed in order to understand the generalization properties as |V | increases. Nevertheless, we believe our work opens up interesting questions, namely better understanding how our results on the energy landscape depend upon specific signal-to-noise ratios, or whether the network parameters can be interpreted mathematically. This could be useful in the study of computational-to-statistical gaps, where our model can be used to inquire about the form of computationally tractable approximations. Another current limitation of our model is that it presumes a fixed number of communities to be detected. Other directions of future research include the extension to the case where the number of communities is unknown and varied, or even increasing with |V |, as well as applications to ranking and edge-cut problems.

Lemma A.1.

|L n (β l ) -L n (β g )| ≤ |L n (β l ) -Ln (β l )| + | Ln (β l ) -Ln ( βg )| + | Ln (β g ) -L n (β g )| (7)

Proof of Lemma A.1. We consider two separate cases: The first case is when Ln (β l ) ≥ Ln (β g ).

Then Ln (β l ) -Ln ( βg ) ≥ Ln (β l ) -Ln (β g ) ≥ 0, and so

|L n (β l ) -L n (β g )| ≤ |L n (β l ) -Ln (β l )| + | Ln (β l ) -Ln ( βg )| + | Ln (β g ) -L n (β g )|.

The other case is when Ln (β l ) < Ln (β g ). Note that

L n (β l ) ≥ L n (β g ). Then |L n (β l ) -L n (β g )| ≤ |L n (β l ) -Ln (β l )| + | Ln (β g ) -L n (β g )| ≤ |L n (β l ) -Ln (β l )| + | Ln (β l ) -Ln ( βg )| + | Ln (β g ) - L n (β g )|.

Hence, to bound the "energy gap" |L n (β l ) -L n (β g )|, if suffices to bound the three terms on the right hand side of Lemma A.1 separately. First, we consider the second term, | Ln (β l ) -Ln ( βg )|.

Let

γ l = R T n β l , γ g = R T n β g and γg = R T n βg . Define S n (γ) = L n (R -T n γ) and Sn (γ) = Ln (R -T

n γ), for any γ ∈ R M . Thus, we apply a change-of-variable and try to bound

EQUATION

Next, we relate the left hand side of the inequality above to cos(γ l , γg ), thereby obtaining an upper bound on [1 -cos 2 (γ l , γg )], which will then be used to bound | Sn (γ l ) -Sn (γ g )|.

EQUATION

Thus, if we define

EQUATION

To bound λ 1 (∇ 2 Sn (γ)), we bound λ 1 (Q 1 ) and Q 2 as follows:

Since Ān is symmetric, let γ1 , . . . γM be the orthonormal eigenvectors of Ān corresponding to nonincreasing eigenvalues l 1 , . . . l M . Note that the global minimum satisfies γg = ±γ 1 . Write

γ = M i=1 α i γi , and let ᾱi = αi √ M i=1 α 2 i . Then | cos(γ, γg )| = | cos(γ, γ1 )| = |ᾱ 1 |.

Then,

EQUATION

To bound Q 2 :

EQUATION

Thus,

EQUATION

This yields the desired lemma.

Combining inequality 8 and Lemma A.3, we get

EQUATION

Thus, to bound the angle between γ l and γg , we can aim to bound ∇ Sn (γ l ) and ∇ 2 S n (γ l ) -∇ 2 Sn (γ l ) as functions of the quantities µ n , ν n and δ n .

Lemma A.4.

EQUATION

Proof of Lemma A.4.

EQUATION

Combining equations 17 and 18, we get

EQUATION

Then, by the generalized Hölder's inequality,

∇ Sn (γ l ) ≤ 2 γ l E|λ 1 (A n )| 3 E ∆B n 3 E 1 |λ M (B n )| 3 1 3 + 3 E|λ 1 (A n )| 3 E ∆B n 3 E 1 |λ M (B n )| 6 1 3 + E|λ 1 (A n )| 3 E ∆B n 6 E 1 |λ M (B n )| 6 1 3 . (21)

Hence, written in terms of the quantities µ n , ν n and δ n , we have

EQUATION

and for k ∈ {1, 2, 3, 4, 5}, if X is a nonnegative random variable,

EX k =1 X>1 EX k + 1 X≤1 EX k ≤1 + 1 X≤1 EX 6 ≤1 + EX 6 (25) Therefore, E|λ 1 (B n )| 6 ≤ 64 + 63E ∆B n 6 .

From now on, for simplicity, we introduce

EQUATION

Proof of Lemma A.6.

EQUATION

where

EQUATION

EQUATION

EQUATION

H 4 = (γ T γ) 2 A n γγ T B n -(γ T B n γ) 2 Aγγ T (γ T B n γ) 2 (γ T γ) 2 (31) Thus, ∇ 2 S n (γ) -∇ 2 Sn (γ) ≤ 2E H 1 + 2E H 2 + 8E H 3 + 8E H 4

, and we try to bound each term on the right hand side separately.

For the first term, there is

EQUATION

Applying generalized Hölder's inequality, we obtain

EQUATION

For the second term, there is

EQUATION

Hence,

EQUATION

Applying generalized Hölder's inequality, we obtain

EQUATION

For H 3 , note that

EQUATION

Hence,

EQUATION

Thus,

EQUATION

Applying generalized Hölder's inequality, we obtain

EQUATION

For the last term,

EQUATION

Thus,

EQUATION

Applying generalized Hölder's inequality, we obtain

EQUATION

Therefore, summing up the bounds above, we obtain

EQUATION

Hence, combining inequality 15, Lemma A.4 and Lemma A.6, we get

EQUATION

Hence, we have proved the theorem, with ηn,µn,νn,δn = 2ηnµnνnδn[2+C(δn,νn)] 1-ηnµnνnδn .

B BACKGROUND B.1 GRAPH MIN-CUTS AND SPECTRAL CLUSTERING:

We consider graphs G = (V, E), modeling a system of N = |V | elements presumed to exhibit some form of community structure. The adjacency matrix A associated with G is the N × N binary matrix such that A i,j = 1 when (i, j) ∈ E and 0 otherwise. We assume for simplicity that the graphs are undirected, therefore having symmetric adjacency matrices. The community structure is encoded in a discrete label vector s : V → {1, . . . , C} that assigns a community label to each node, and the goal is to estimate s from observing the adjacency matrix.

In the binary case, we can set s(i) = ±1 without loss of generality. Furthermore, we assume that the communities are associative, which means two nodes from the same community are more likely to be connected than two nodes from the opposite communities. The quantity i,j

(1 -s(i)s(j))A i,j

measures the cost associated with cutting the graph between the two communities encoded by s, and we wish to minimize it under appropriate constraints (Newman, 2006) . Note that i,j A i,j = s T Ds, with D = diag(A1) (called the degree matrix), and so the cut cost can be expressed as a positive semidefinite quadratic form min

s(i)=±1

s T (D -A)s = s T ∆s that we wish to minimize. This shows a fundamental connection between the community structure and the spectrum of the graph Laplacian ∆ = D -A, which provides a powerful and stable relaxation of the discrete combinatorial optimization problem of estimating the community labels for each node. The eigenvector of ∆ associated with the smallest eigenvalue is, trivially, 1, but its Fiedler vector (the eigenvector associated with the second smallest eigenvalue) reveals important community information of the graph under appropriate conditions (Newman, 2006) , and is associated with the graph conductance under certain normalization schemes (Spielman, 2015) .

Given linear operator L(A) extracted from the graph (that we assume symmetric), we are thus interested in extracting eigenvectors at the edge of its spectrum. A particularly simple algorithm is the power iteration method. Indeed, the Fiedler vector of L(A) can be obtained by first extracting the leading eigenvector v of Ã = L(A) I -L(A), and then iteratively compute

y (n) = Ãw (n-1) , w (n) = y (n) -y (n) , v v y (n) -y (n) , v v .

Unrolling power iterations and recasting the resulting model as a trainable neural network is akin to the LISTA sparse coding model, which unrolled iterative proximal splitting algorithms (Gregor & LeCun, 2010) .

Despite the appeal of graph Laplacian spectral approaches, it is known that these methods fail in sparsely connected graphs (Krzakala et al., 2013) . Indeed, in such scenarios, the eigenvectors of the graph Laplacian concentrate on nodes with dominant degrees, losing their correlation with the community structure. In order to overcome this important limitation, people have resorted to ideas inspired from statistical physics, as explained next.

B.2 PROBABILISTIC GRAPHICAL MODELS AND BELIEF-PROPAGATION:

Graphs with labels on nodes and edges can be cast as a graphical model where the aim of clustering is to optimize label agreement. This can be seen as a posterior inference task. If we simply assume the graphical model is a Markov Random Field (MRF) with trivial compatibility functions for cliques greater than 2, the probability of a label configuration σ is given by

EQUATION

Generally, computing marginals of multivariate discrete distributions is exponentially hard. For instance, in the case of P(σ i ) we are summing over |X| n-1 terms (where X is the state space of discrete variables). But if the graph is a tree, we can factorize the MRF more efficiently to compute the marginals in linear time via a dynamic programming method called the sum-product algorithm, also known as belief propagation (BP). An iteration of BP is given by

EQUATION

The beliefs (b i→j (σ i )) are interpreted as the marginal distributions of σ i . Fixed points of BP can be used to recover marginals of the MRF above. In the case of the tree, the correspondence is exact:

P i (σ i ) = b i (σ i ).

Certain sparse graphs, like SBM with constant degree, are locally similar to trees for such an approximation to be successful (Mossel et al., 2014) . However, convergence is not guaranteed in graphs that are not trees. Furthermore, in order to apply BP, we need a generative model and the correct parameters of the model. If unknown, the parameters can be derived using expectation maximization, further adding complexity and instability to the method since it is possible to learn parameters for which BP does not converge.

B.3 NON-BACKTRACKING OPERATOR AND BETHE HESSIAN:

The BP equations have a trivial fixed-point where every node takes equal probability in each group. Linearizing the BP equation around this point is equivalent to spectral clustering using the nonbacktracking matrix (NB), a matrix defined on the directed edges of the graph that indicates whether two edges are adjacent and do not coincide. Spectral clustering using NB gives significant improvements over spectral clustering with different versions of the Laplacian matrix L and the adjacency matrix A. High degree fluctuations drown out the signal of the informative eigenvalues in the case of A and L, whereas the eigenvalues of NB are confined to a disk in the complex plane except for the eigenvalues that correspond to the eigenvectors that are correlated with the community structure, which are therefore distinguishable from the rest.

NB matrices are still not optimal in that they are matrices on the edge set and also asymmetric, therefore unable to enjoy tools of numerical linear algebra for symmetric matrices. Saade et al. (2014) showed that a spectral method can do as well as BP in the sparse SBM using the Bethe Hessian matrix defined by BH(r) := (r 2 -1)I -rA + D, where r is a scalar parameter. This is due to a one-to-one correspondence between the fixed points of BP and the stationary points of the Bethe free energy (corresponding Gibbs energy of the Bethe approximation) (Saade et al., 2014) . The Bethe Hessian is a scaling of the Hessian of the Bethe free energy at an extrema corresponding to the trivial fixed point of BP. Negative eigenvalues of BH(r) correspond to phase transitions in the Ising model where new clusters become identifiable. The success of the spectral method using the Bethe Hessian gives a theoretical motivation for having a family of matrices including I, D and A in our GNN defined in Section 4, because in this way the GNN is capable of expressing the algorithm of performing power iteration on the Bethe Hessian. While belief propagation requires a generative model, and the spectral method using the Bethe Hessian requires the selection of the parameter r, whose optimal value also depends on the underlying generative model, the GNN does not need a generative model and is able to learn and then make predictions in a data-driven fashion. motivates other estimation algorithms based on motif-counting that require knowledge of the model likelihood function (Sankararaman & Baccelli, 2018) .

Table 3 shows the performance of GNN and LGNN on the binary GBM model, obtained with d = 2, n = 500, T = 5 √ 2 and varying S, as well as the performances of two spectral methods, using respectively the normalized Laplacian and the Bethe Hessian, which approximates BP around its stationary solution. We note that LGNN model, thanks to its added flexibility and the multiscale nature of its generators, is able to significantly outperform both spectral methods as well as the baseline GNN.

C.2 MIXTURE OF BINARY SBM:

We report here our experiments on the SBM mixture, generated with G ∼ SBM (n = 1000, p = k d -q, q ∼ Unif(0, dd), C = 2) , where the average degree d is either fixed constant or also randomized with d ∼ Unif(1, t). Figure 4 shows the overlap obtained by our model compared with several baselines. Our GNN model is either competitive with BH or outperforms BH, which achieves the state of the art along with BP Saade et al. (2014) , despite not having any access to the underlying generative model (especially in cases where GNN was trained on a mixture of SBM and thus must be able to generalize the r parameter in BH). They all outperform by a wide margin spectral clustering methods using the symmetric Laplacian and power method applied to BH I -BH using the same number of layers as our model. Thus GNN's ability to predict labels goes beyond approximating spectral decomposition via learning the optimal r for BH(r). The model architecture could allow it to learn a higher dimensional function of the optimal perturbation of the multiscale adjacency basis, as well as nonlinear power iterations, that amplify the informative signals in the spectrum.

C.3 ABLATION STUDIES OF GNN AND LGNN:

Compared to f , each of h, i and k has one fewer operator in F, and j has two fewer. We see that with the absence of A (2) , k has much worse performance than the other four, indicating the importance of the power graph adjacency matrices. Interestingly, with the absence of I, i actually has better average accuracy than f . One possibly explanation is that in SBM, each node has the same expected degree, and hence I may be not very far from D, which might make having both I and D in the family redundant to some extent.

Comparing GNN models a, b and c, we see it is not the case that having larger J will always lead to better performance. Compared to f , GNN models c, d and e have similar numbers of parameters but all achieve worse average test accuracy, indicating that the line graph structure is essential for the good performance of LGNN in this experiment. In addition, l also performs worse than f , indicating the significance of the non-backtracking line graph compared to the symmetric line graph.

:

In this work, however, we assume that C is fixed.

The term "spatial batch normalization" may be slightly misleading, since each batch only contain one graph in our experiments. However, we used spatial batch normalization in our implementation for convenience, as it performs orthogonalization and normalization, as further explained in the next paragraph.

to account for the invariance up to global flip of label

Code is available at https://github.com/zhengdao-chen/GNN4CD

Our implementation of GAT is based on https://github.com/Diego999/pyGAT. We modified the code so that the number of layers in the network is flexible, and also added spatial batch normalization at the end of each layer, similar to in our GNN and LGNN, as our experiments showed that including spatial batch normalization improved the performance. The results in sections

6.1 and 6.2 are from a trained GAT with 30 layers and 8 features.

Figures and Tables:

Figure 2. Construction of the line graph L(G) using the non-backtracking Operator. The nodes of L(G) correspond to oriented edges of G.

and and finally use (θ) = inf π∈S C -i∈V log ōi,π( ȳi) as an approximate loss value, which only involves a permutation group of size ( C!).

Figure 3. Binary associative SBM detection (C = 2, p > q). X-axis corresponds to SNR, and Y-axis to overlap between the prediction and the ground truth. GNN LGNN LGNN-L LGNN-S GAT BP Avg. 0.18 0.21 0.18 0.18 0.16 0.14 Std. Dev. 0.04 0.05 0.04 0.04 0.04 0.02

Figure 4. GNN mixture (Graph Neural Network trained on a mixture of SBM with average degree 3), GNN full mixture (GNN trained over different SNR regimes), BH( √ d) and BH(-√ d). left: k = 2. We verify that BH(r) models cannot perform detection at both ends of the spectrum simultaneously.

Comparison of the node classification accuracy by different models on the three SNAP datasets. Note that the average accuracy was computed graph-wise with each graph weighted by its size, while the standard deviation was computed graph-wise with equal weights among the graphs.

Table content: <table><tr><td/><td>train/test</td><td colspan="2">Avg |V| Avg |E|</td><td/><td colspan="5">GNN LGNN LGNN-S GAT AGMfit</td></tr><tr><td colspan="2">Amazon 805/142</td><td>60</td><td>161</td><td>Avg.</td><td>0.97</td><td>0.96</td><td>0.97</td><td>0.95</td><td>0.90</td></tr><tr><td/><td/><td/><td/><td colspan="2">Std. Dev. 0.12</td><td>0.13</td><td>0.11</td><td>0.13</td><td>0.13</td></tr><tr><td>DBLP</td><td>4163/675</td><td>26</td><td>77</td><td>Avg.</td><td>0.90</td><td>0.90</td><td>0.89</td><td>0.88</td><td>0.79</td></tr><tr><td/><td/><td/><td/><td colspan="2">Std. Dev. 0.13</td><td>0.13</td><td>0.13</td><td>0.13</td><td>0.18</td></tr><tr><td colspan="3">Youtube 20000/1242 93</td><td>201</td><td>Avg.</td><td>0.91</td><td>0.92</td><td>0.91</td><td>0.90</td><td>0.59</td></tr><tr><td/><td/><td/><td/><td colspan="2">Std. Dev. 0.11</td><td>0.11</td><td>0.11</td><td>0.13</td><td>0.16</td></tr></table>

The effects of different architectures and choices of the operator family for GNN and LGNN, as demonstrated by their performance on the 5-class dissociative SBM experiments with the exact setup as in Section 6.2. For LGNN, FB is the same as FA except for changing A to B.

Table content: <table><tr><td/><td/><td colspan="3">#layers #features J</td><td>FA</td><td colspan="2">#parameters Avg.</td><td>Std. Dev.</td></tr><tr><td colspan="2">(a) GNN</td><td>30</td><td>8</td><td>2</td><td>I, D, A, A 2</td><td>8621</td><td colspan="2">0.1792 0.0385</td></tr><tr><td colspan="2">(b) GNN</td><td>30</td><td>8</td><td>4</td><td>I, D, A, ..., A (4)</td><td>12557</td><td colspan="2">0.1855 0.0438</td></tr><tr><td colspan="2">(c) GNN</td><td>30</td><td>8</td><td colspan="3">11 I, D, A, ..., A (11) 26333</td><td colspan="2">0.1794 0.0359</td></tr><tr><td colspan="2">(d) GNN</td><td>30</td><td>15</td><td>2</td><td>I, D, A, A (2)</td><td>28760</td><td colspan="2">0.1894 0.0388</td></tr><tr><td colspan="2">(e) GNN</td><td>30</td><td>12</td><td>4</td><td>I, D, A, ..., A (4)</td><td>27273</td><td colspan="2">0.1765 0.0371</td></tr><tr><td colspan="2">(f) LGNN</td><td>30</td><td>8</td><td>2</td><td>I, D, A, A (2)</td><td>25482</td><td colspan="2">0.2073 0.0481</td></tr><tr><td colspan="3">(g) LGNN-L 30</td><td>8</td><td>2</td><td>I, D, A, A (2)</td><td>25482</td><td colspan="2">0.1822 0.0395</td></tr><tr><td colspan="2">(h) LGNN</td><td>30</td><td>8</td><td>2</td><td>I, A, A (2)</td><td>21502</td><td colspan="2">0.1981 0.0529</td></tr><tr><td>(i)</td><td>LGNN</td><td>30</td><td>8</td><td>2</td><td>D, A, A (2)</td><td>21502</td><td colspan="2">0.2212 0.0581</td></tr><tr><td>(j)</td><td>LGNN</td><td>30</td><td>8</td><td>2</td><td>A, A (2)</td><td>17622</td><td colspan="2">0.1954 0.0441</td></tr><tr><td colspan="2">(k) LGNN</td><td>30</td><td>8</td><td>1</td><td>I, D, A</td><td>21502</td><td colspan="2">0.1673 0.0437</td></tr><tr><td>(l)</td><td colspan="2">LGNN-S 30</td><td>8</td><td>2</td><td>I, D, A, A (2)</td><td>21530</td><td colspan="2">0.1776 0.0398</td></tr></table>

Peer review (iclr2019):
 This paper presents a study of the community detection problem via graph neural networks. The presented results open the possibility that neural networks are able to discover the optimal algorithm for a given task. This is rather convincingly demonstrated on the example of the stochastic block model, where the optimal performance is known (for 2 symmetric groups) or strongly conjectured (for more groups). The method is rather computationally demanding, and also somewhat unrealistic in the aspect that the training examples might not be available, but for a pioneering study of this kind this is well acceptable. Despite my overall very positive opinion, I found a couple of claims that are misleading and overall hurt the quality of the paper, and I would strongly suggest to the authors to adjust these claims: ** The method is claimed to "even improve upon current computational thresholds in hard regimes." This is misleading, because (as correctly stated in the body of the paper) the computational threshold to which the paper refers apply in the limit of large graph sizes whereas the observed improvements are for finite sizes. It is shown here that for finite sizes the present method is better than belief propagation. But this clearly does not imply that it improves the conjectured computational thresholds that are asymptotic. At best this is an interesting hypothesis for future work, not more. ** The energy landscape is analyzed "under certain simplifications and assumptions". Conclusions state "an interesting transition from rugged to simple as the size of the graphs increase under appropriate concentration conditions." This is very vague. It would be great if the paper could offer intuitive explanation of there simplifications and assumptions that is between these unclear remarks and the full statement of the theorem and the proof that I did not find simple to understand. For instance state the intuition on in which region of parameters are those results true and in which they are not. ** "multilinear fully connected neural networks whose landscape is well understood (Kawaguchi, 2016)." this is in my opinion grossly overstated. While surely that paper presents interesting results, they are set in a regime that lets a lot to be still understood about landscape of fully connected neural networks. It is restricted to specific activation functions, and the results for non-linear networks rely on unjustified simplifications, the sample complexity trade-off is not considered, etc. Misprint: Page 2: cetain -> certain.