\section{Introduction}\label{sec:introduction}

In recent years there has been a surge of research interest in developing varieties of Graph Neural Networks (GNNs) -- specialized deep learning architectures for dealing with graph-structured data, such as social networks \cite{perozzideepwalk}, recommender graphs \cite{ying2018graph}, or molecular graphs~\cite{defferrard2016,ying2018}.
GNNs leverage the structure of the data as computational graph, allowing the information to propagate across the edges of graphs~\cite{scarselli2008}.
When many real-wold systems are represented as graphs, they exhibit locally inhomogeneous distributions of edges, forming \emph{clusters} (also called \emph{communities} or \emph{modules}) -- groups of nodes with high in-group edge density, and relatively low out-group density. Clusters can correspond to interesting phenomena in the underlying graph, for example to education \cite{traud2011comparing} or employment \cite{papacharissi2009virtual} in social graphs.
GNNs have been shown to benefit from leveraging higher-order structural information that could arise from clusters
~\cite{chen2018harp,ying2018,lee2019}, for example through pooling or trainable attention over edges~\cite{velickovic2017}.

Interestingly, most existing work on GNNs to leverage higher-order structure does not directly address node partitioning or the estimation of clusters within the computational graph. Furthermore, most works explore these mechanisms only within a semi-supervised or supervised framework, ignoring the fact that \emph{unsupervised} graph clustering is often an extremely useful end-goal in itself -- whether for data exploration \cite{perozzi2018discovering},  visualization \cite{clemenccon2012hierarchical, cui2008geometry}, genomic feature discovery \cite{cabreros2016detecting}, anomaly detection \cite{perozzi2016scalable}, or for many other use-cases discussed e.g.\ in \cite{fortunato2016community}. Additionally, many of the existing structure-aware methods have undesirable properties, such as relying on a multi-step optimization process which does not allow to optimize the objective via gradient descent end-to-end \cite{perozzi2014focused}.

\clearpage
In this work, we take an \emph{ab initio} approach to the clustering problem in the GNN domain, bridging the gap between traditional graph clustering objectives and deep neural networks.
We start by drawing a connection between graph pooling, which was typically studied in the literature as a regularizer for supervised GNN architectures, and fully unsupervised clustering.
Specifically, we contribute:
\vspace{-2mm}
\begin{itemize}[leftmargin=0.2cm,itemindent=.3cm,labelwidth=\itemindent,labelsep=0cm,align=left,topsep=0pt,itemsep=-1ex,partopsep=0ex,parsep=1ex]
    \item \textsc{DMoN}\xspace, a method for unsupervised clustering in GNNs that allows optimization of cluster assignments in an end-to-end differentiable way.
    \item An empirical study of performance on synthetic graphs, illustrating the problems with existing work and how \textsc{DMoN}\xspace allows for improved model performance in those regimes.
    \item Thorough experimental evaluation on real-world data, showing that many pooling methods poorly reflect hierarchical structures and are not able to make use of both graph structure and node attributes.
\end{itemize}
 %
\section{Related Work}\label{sec:related-work}

\begin{table}[!t]
\small
\centering{
\newcolumntype{R}{>{\raggedleft\arraybackslash}X}
\newcolumntype{C}{>{\centering\arraybackslash}X}
\newcolumntype{L}{>{\hsize=1.05\hsize}C}
\newcolumntype{S}{>{\hsize=0.75\hsize}C}
\caption{Related work in terms of six desirable clustering properties outlined in Section~\ref{sec:related-work}.}
\label{tbl:related-work}
\begin{tabularx}{\linewidth}{p{1cm}SSSLLS@{}}
\toprule
\emph{method} & Trainable & Unsupervised & Sparse & \mbox{Node pooling} & \mbox{Soft assignments} & Stable \\
\midrule
Graclus~\cite{defferrard2016} & \textcolor{cycle2}{\ding{56}} & \cmark & \cmark & \cmark & \textcolor{cycle2}{\ding{56}} & \cmark \\
DiffPool~\cite{ying2018} & \cmark & \cmark & \textcolor{cycle2}{\ding{56}} & \cmark & \cmark & \textcolor{cycle2}{\ding{56}} \\
\mbox{Top-k}~\cite{gao2019} & \cmark & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \cmark \\
\mbox{SAG}~\cite{lee2019} & \cmark & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \textcolor{cycle2}{\ding{56}} & \cmark \\
\mbox{MinCut}~\cite{bianchi2020} & \cmark & \cmark & \cmark & \cmark & \cmark & \textcolor{cycle2}{\ding{56}} \\
\midrule
\mbox{\textsc{DMoN}\xspace} & \cmark & \cmark & \cmark & \cmark & \cmark & \cmark \\
\bottomrule
\end{tabularx}}
\end{table} \vspace{-4mm}

Our work builds upon a rich line of research on graph neural networks and graph pooling methods.

\spara{Graph Neural Networks} (GNNs)~\cite{scarselli2008,duvenaud2015,niepert2016,gilmer2017,kipf2017} allow end-to-end differentable losses over data with arbitrary structure.  They have been applied to an incredible range of applications, from social networks \cite{perozzideepwalk}, to recommender systems \cite{ying2018graph}, to computational chemistry~\cite{gilmer2017}.
While GNNs are flexible enough to allow for unsupervised losses, most work follows the semi-supervised setting for node classification from \cite{kipf2017}.
For a complete introductions to the vast topic we refer interested readers to detailed surveys~\cite{bronstein2017,hamilton2017survey,chami2020}.

\spara{Unsupervised training of GNNs} is commonly done via maximizing mutual information~\cite{belghazi2018,hjelm2018,tschannen2019,song2019} in a self-supervised fashion.
Deep Graph Infomax (DGI)~\cite{velickovic2018} adapted the mutual informaton-based learning from Deep InfoMax~\cite{hjelm2018}, learning unsupervised representations for nodes in attributed graphs.
InfoGraph~\cite{sun2020} extended the idea to learning representations of whole graphs instead of just nodes. A very similar approach was introduced independently in~\cite{hu2020} in the context of pre-training GNNs for producing representations of graphs, which was first tackled in~\cite{navarin2018}.

\spara{Graph pooling} aims to tackle the hierarchical nature of graphs via iterative coarsening.
Early architectures~\cite{defferrard2016} resorted to fixed axiomatic pooling, with no optimization of clustering while the network learns.
DiffPool~\cite{ying2018} suggests to include a \emph{learnable} pooling to GNN architecture. To help the convergence, DiffPool includes a link prediction loss to help encapsulate the clustering structure of graphs and an additional entropy loss %
to penalize soft assignments. Top-k~\cite{gao2019} and SAG pooling~\cite{lee2019} learn to sparsify the graph (select top-k edges for each node) with learned weights. MinCutPool~\cite{bianchi2020} pooling studies differentiable formulation of spectral clustering as a pooling strategy.

We summarize mainstream graph pooling methods in Table~\ref{tbl:related-work} in terms of six desirable properties related to their clustering capabilities:
\vspace{-2mm}
\begin{itemize}[leftmargin=0.2cm,itemindent=.3cm,labelwidth=\itemindent,labelsep=0cm,align=left,topsep=0pt,itemsep=-1ex,partopsep=0ex,parsep=1ex]
    \item \textbf{Trainable.} End-to-end training allows to capture both graph structure and node features.
    \item \textbf{Unsupervised} training is a desirable setting for clustering models. Works on \emph{supervised} graph clustering~\cite{yang2019,wang2019cvpr} are outside of the scope of this work.
    \item \textbf{Sparse.} As graphs in the real-world vary in size and sparsity, methods can not be limited by a $\mathcal{O}(n^2)$ link prediction objectives, like DiffPool~\cite{ying2018}, or computing $\mathbf{A}^2$, like top-k pooling methods~\cite{gao2019,lee2019} .
    \item \textbf{Node aggregation} is crucial for our interpretation of graph pooling in terms of graph clustering. Both Top-k~\cite{gao2019} and SAG pooling~\cite{lee2019} only sparsify the graph and do not reduce the nodeset.
    \item \textbf{Soft assignments} allow for more flexibility reasoning about the interactions of clusters.
    \item \textbf{Stable} -- the method should be stable in terms of the graph structure. DiffPool~\cite{ying2018} is not stable in terms of graph sparsity, while MinCutPool~\cite{bianchi2020} can not deal with uneven degree distribution.
\end{itemize}

\spara{Graph embeddings}~\cite{perozzi2014,grover2016,tsitsulin2018} can be thought of as (very restricted) unsupervised GNNs with an identity feature matrix, meaning each node learns its own positional representation~\cite{you2019}. The learning process in graph embeddings is often done in a similar way to DGI though noise contrastive estimation~\cite{gutmann2010}.
As far as we know, all pooling strategies for learning node embeddings without attributes have been axiomatic~\cite{chen2018,liang2018,deng2020}. %
\section{Preliminaries}\label{sec:preliminaries}

We introduce the necessary background for \textsc{DMoN}\xspace{}, starting with the problem formulation, reviewing common graph clustering objectives and how they can be made differentiable efficiently.

Graph $G=(V,E)$ is defined via a set of nodes $V=(v_1,\ldots,v_n), |V|=n$ and edges $ E\subseteq V\times V, |E|=m$.
We are interested in measuring the quality of graph partitioning function $\mathcal{F}: n \rightarrow k$ that splits the set of nodes $V$ into $k$ partitions $V_i=\{v_j, \mathcal{F}(v_j)=i\}$.
Additionally, in contrast to standard graph clustering, we are also provided with node attributes $\mathbf{X}\in\mathbb{R}^{n\times s}$ that provide additional information not reflected in the graph structure, but also correlated with it.

\subsection{Graph Clustering Quality Functions}\label{ssec:modularity}

Design of the objective function is crucial for the algorithm performance.
We review two families of clustering quality functions amenable to spectral optimization, and review some of their shortcomings.

\spara{Cut-based metrics.} In his seminal work~\cite{fiedler1973algebraic}, Fiedler suggested that the second (Fiedler) eigenvector of a graph Laplacian produces a graph \emph{cut} minimal in terms of the weight of the edges.
This plain notion of cut degenerates on real-world graphs, as it does not require partitions to be balanced in terms of size.
It is possible to get normalized partitions with the use of ratio cut~\cite{wei1989towards}, which normalizes the cut by the product of the number of nodes in two partitions, or normalized cut~\cite{shi2000normalized}, which uses total edge volume of the partition as normalization.

In real networks, however, there is evidence \emph{against} existence of good cuts~\cite{leskovec2008statistical} in ground-truth communities.
This can be explained by the fact than a single node implicitly participates in many different clusters~\cite{epasto2017ego}, e.g.\ a person in a social network is simultaneously connected with family and work friends, forcing the algorithm to merge these communities together.

Recently, MinCutPool~\cite{bianchi2020} adapted the notion of the normalized cut to use as a regularizer for pooling.
While, theoretically, MinCutPool objective should be suitable for clustering nodes in graphs, we (i) show that it does not optimize its own objective function and (ii) provide evidence against this in the synthetic and real-world experiments.

\spara{Modularity}~\cite{newman2006} approaches the same problem from a statistical perspective, incorporating a \emph{null model} to quantify the significance of the clustering with respect to the random graph.
In a fully random graph with given degrees, nodes $u$ and $v$ with degrees $d_u$ and $d_v$ are connected with probability $\nicefrac{d_ud_v}{2m}$.
Modularity measures the divergence between the intra-cluster edges from the expected one:
\begin{equation}\label{eq:modularity}
    \mathcal{Q} = \frac{1}{2m}\sum_{ij}\left[\mathbf{A}_{ij}-\frac{d_id_j}{2m}\right]\delta(c_i, c_j),
\end{equation}
\noindent where $\delta(c_i, c_j)=1$ if $i$ and $j$ are in the same cluster and 0 otherwise.
Note $\mathcal{Q}\in(\nicefrac{-1}{2}; 1]$ (it is 0 when there is no correlation of clusters with edge density), but it is not necessarily maximized at 1, and is only comparable across graphs with the same degree distribution.
While problems with the modularity metric have been identified \cite{good2010performance}, it remains one of the most commonly-used and eminently useful graph clustering metrics in scientific literature \cite{fortunato2016community}.

\subsection{Spectral Modularity Maximization}\label{ssec:spectral-modularity}
Maximizing the modularity is proven to be $\mathbf{NP}$-hard\xspace~\cite{brandes2006maximizing}, however, a spectral relaxation of the problem can be solved efficiently~\cite{newman2006finding}.
Let $\mathbf{C}\in{0, 1}^{n\times k}$ be the cluster assignment matrix and $\mathbf{d}$ be the degree vector.
Then, with \emph{modularity matrix} $\mathbf{B}$ defined as $\mathbf{B}=\mathbf{A}-\frac{\mathbf{d}\vD^\top}{2m}$, the modularity $\mathcal{Q}$ can be reformulated as:
\begin{equation}\label{eq:spectral_modularity}
    \mathcal{Q} = \frac{1}{2m}\tr(\mathbf{C}^\top\mathbf{B}\mathbf{C})
\end{equation}
Relaxing $\mathbf{C}\in\mathbb{R}^{n\times k}$, the optimal $\mathbf{C}$ maximizing $\mathcal{Q}$ is the top-k eigenvectors of the modularity matrix $\mathbf{B}$.
While $\mathbf{B}$ is dense, iterative eigenvalue solvers can take advantage of the fact that $\mathbf{B}$ is a sum of a sparse $\mathbf{A}$ and rank-one matrix $-\frac{\mathbf{d}\vD^\top}{2m}$, meaning that the matrix-vector product $\mathbf{B}\mathbf{x}$ can be computed efficiently as
\begin{equation*}\label{eq:modularity-matvec}
    \mathbf{B}\mathbf{x} = \mathbf{A}\mathbf{x} - \frac{\mathbf{d}^\top\mathbf{x}\mathbf{d}}{2m}
\end{equation*}

\noindent and optimized efficiently with iterative methods such as power iteration or Lanczos algorithm.
One can then obtain clusters by means of spectral bisection~\cite{newman2006finding} with iterative refinement akin to Kernighan-Lin algorithm~\cite{kernighan1970efficient}.
However, these formulations operate entirely on the graph structure, and it is non-trivial to adapt them to work with attributed graphs.

\subsection{Graph Neural Networks}\label{sec:gnns}

Graph Neural Networks are a flexible class of models that perform nonlinear feature aggregation with respect to graph structure.
For the purposes of this work, we consider transductive GNNs that output a single embedding per node.
Graph convolutional networks (GCNs)~\cite{kipf2017} are simple yet effective~\cite{shchur2018} message-passing networks that fit our criteria.
Let $\mathbf{X}^0\in\mathbb{R}^{n\times s}$ be the initial node features and $\Tilde{\mathbf{A}} = \mathbf{D}^{-\frac{1}{2}}\mathbf{A}\mathbf{D}^{-\frac{1}{2}}$ be the normalized adjacency matrix, the output of $t$-th layer $\mathbf{X}^{t+1}$ is
\begin{equation}\label{eq:gcn}
    \mathbf{X}^{t+1} = \mathrm{SeLU}(\Tilde{\mathbf{A}}\mathbf{X}^{t}\mathbf{W} + \mathbf{X}\mathbf{W}_{\mathrm{skip}})
\end{equation}

We make two changes to the classic GCN architecture: first, we remove the self-loop creation and instead use an $\mathbf{W}_{\mathrm{skip}}\in\mathbb{R}^{s\times s}$ trainable skip connection, and, second, we replace ReLU nonlinearity with SeLU~\cite{klambauer2017self} for better convergence. %
\section{Method}\label{sec:method}

In this section, we present \textsc{DMoN}\xspace, our method for attributed graph clustering with graph neural networks.
Inspired by the modularity function and its spectral optimization, we propose a fully differentiable unsupervised clustering objective which optimizes soft cluster assignments using a null model to control for inhomogeneities in the graph.
We then discuss the challenge of regularizing cluster assignments, and present collapse regularization, a softer version of orthogonality loss, that against trivial solutions.

\subsection{\textsc{DMoN}\xspace: Deep Modularity Networks}
The challenge of clustering boils down to defining an optimization procedure over the cluster assignment matrix $\mathbf{C}$.
In \textsc{DMoN}\xspace, we propose to obtain $\mathbf{C}$ via the output of a softmax function, which allows the (soft) cluster assignment to be differentiable.
The input to the cluster assignment can be any differentiable message passing function, but here we specifically consider the case where a graph convolutional network is used to obtain soft clusters for each node as follows:
\begin{equation}\label{eq:dmon-architecture}
    \mathbf{C} = \mathrm{softmax}(\mathrm{GCN}(\Tilde{\mathbf{A}}, \mathbf{X})),
\end{equation}

\noindent where $\mathrm{GCN}$ is a (possibly) multi-layer convolutional network operating on an normalized adjacency matrix $\Tilde{\mathbf{A}} = \mathbf{D}^{-\frac{1}{2}}(\mathbf{A})\mathbf{D}^{-\frac{1}{2}}$.

We then propose to optimize this assignment with the following objective, which combines insights from spectral modularity maximization \eqref{eq:spectral_modularity} with a regularization to ensure informative clusters:

\begin{equation}\label{eq:dmon-objective}
    \mathcal{L}_\text{{\textsc{DMoN}\xspace}} = \underbrace{-\frac{1}{2m}\tr ( \mathbf{C}^\top \mathbf{B} \mathbf{C})}_{\mathrm{modularity}} +
    \underbrace{\frac{\sqrt{k}}{n}\left\lVert \sum_i\mathbf{C}^\top_i\right\rVert_F - 1}_{\mathrm{collapse~regularization}},
\end{equation}

where $\lVert\cdot\rVert_F$ is the Frobenius norm.
We decompose the computation of $\tr(\mathbf{C}^\top\mathbf{B}\mathbf{C})$ as a sum of sparse matrix-matrix multiplication and rank-one degree normalization $\tr(\mathbf{C}^\top\mathbf{A}\mathbf{C}-\mathbf{C}^\top\mathbf{d}^\top\mathbf{d}\mathbf{C})$. This allows us to efficiently optimize \textsc{DMoN}\xspace parameters in sparse regime.

\subsection{Collapse regularization}\label{ssec:clustersize}

\begin{wrapfigure}[16]{R}{0.4\textwidth}
\vspace{-12mm}
\centering
\begin{tikzpicture}
\begin{groupplot}[group style={
                      group name=myplot,
                      group size= 1 by 2,vertical sep=0.75cm},
                      height=3cm,
                      width=1.2\linewidth,
                      title style={at={(0.5,0.9)},anchor=south},
                      every axis x label/.style={at={(axis description cs:0.5,-0.2)},anchor=north},
                      yticklabels={,,},
                      xmin=0,
                      xmax=200]
\nextgroupplot[
 	title = \textbf{\textsc{DMoN}\xspace{}},
 	legend columns=2,
	legend style={at={(0.5,1.4)},anchor=south},
    legend entries={Objective, Regularizer},
	xticklabels={,,},
]
\addplot[very thick,color=cycle2] table[x=epoch,y=mgcn] {data/loss_histories.tex};
\addplot[very thick,color=cycle3] table[x=epoch,y=mgcn-reg] {data/loss_histories.tex};
\nextgroupplot[
 	title = \textbf{MinCutPool},
	xlabel=epoch,
]
\draw [color=black, dashed, very thick, draw opacity=0.75] (0, -0.9972) -- (200,-0.9972);
\draw[->,thick,color=black](axis cs:60,-0.2)--(axis cs:75,-0.8);
\node[anchor=east] (text) at (axis cs:60,-0.2){\tiny \textbf{not improving}};
\addplot[very thick,color=cycle2] table[x=epoch,y=mincut] {data/loss_histories.tex};
\addplot[very thick,color=cycle3] table[x=epoch,y=mincut-reg] {data/loss_histories.tex};
\end{groupplot}
\end{tikzpicture}
\vspace*{-7mm}
\caption{\label{fig:loss-evolution}\small Optimization progress of {MinCut} and \textsc{DMoN}\xspace{} on Cora dataset. MinCut optimizes the regularizer, while \textsc{DMoN}\xspace{} minimizes its main objective.}
\end{wrapfigure}
Without additional constraints on the assignment matrix $\mathbf{C}$, spectral clustering for both min-cut and modularity objectives has spurious local minima: assigning all nodes to the same cluster produces a trivial locally optimal solution that traps gradient-based optimization methods.
MinCut pooling~\cite{bianchi2020} addresses this problem by adapting spectral orthogonality constraint in the form of soft-orthogonality~\cite{bansal2018can} regularization $\left\lVert \mathbf{C}^\top\mathbf{C} - \mathbf{I}\right\rVert_F$.
We notice that this term is overly restrictive when combined with softmax class assignment -- intuitively, when the value range of $\mathbf{C}$ is restricted to $\mathbb{R}\cap[0, 1]$ the optimization of the soft-orthogonality regularizer dominates the loss.

We illustrate this problem in Figure~\ref{fig:loss-evolution}, which depicts the progress of optimization for both methods in terms of their main objective and regularizer term over the course of 200 epochs on Cora dataset.
The soft orthogonality regularization term dominates the optimization for MinCutPool, such that the cut objective becomes \emph{worse than random} over the course of training.

We fix this problem by proposing a relaxed notion of \emph{collapse regularization} that prevents the trivial partition while not restricting the optimization of the main objective.
The regularizer is a Frobenius norm of the (soft) cluster membership counts, normalized to range $[0, 1]$.
It gets value of $0$ when cluster sizes are perfectly balanced, and $1$ in the case all clusters collapse to one.
We also propose to stabilize the training by applying dropout~\cite{srivastava2014dropout} to GNN representations before the softmax, preventing the gradient descent from getting stuck in the local optima of the highly non-convex objective function.
 %
\vspace{-5mm}
\section{Experiments}\label{sec:experiments}
\vspace{-4mm}

We assess the clustering performance of \textsc{DMoN}\xspace{} in terms of both graph clustering and label alignment.
We release the implementation of \textsc{DMoN}\xspace{} at this URL\footnote{\href{https://github.com/google-research/google-research/tree/master/graph_embedding/dmon}{\texttt{github.com/google-research/google-research/tree/master/graph\_embedding/dmon}}}.

\begin{wrapfigure}[10]{r}{7cm}
\vspace{-5mm} %
\footnotesize
\centering{
\newcolumntype{C}{>{\centering\arraybackslash}X}
\caption{\label{tbl:datasets}Dataset statistics.}
\vspace*{-2mm}
\begin{tabularx}{7cm}{p{1.5cm}CCCC@{}}
\toprule
\emph{dataset} & Nodes & Edges & Features & Classes \\
\midrule
Cora\xspace & 2708 & 5278 & 1433 & 7 \\
Citeseer\xspace & 3327 & 4614 & 3703 & 6 \\
Pubmed\xspace & 19717 & 44325 & 500 & 3 \\
\mbox{Amazon PC\xspace} & 13752 & 143604 & 767 & 10 \\
\mbox{Amazon Photo\xspace} & 7650 & 71831 & 745 & 8 \\
\mbox{Coauthor CS\xspace} & 18333 & 81894 & 6805 & 15 \\
\mbox{Coauthor PHY\xspace} & 34493 & 247962 & 8415 & 5 \\
\bottomrule
\end{tabularx}}
\end{wrapfigure}

\spara{Datasets.}
We use 7 real-world datasets for assessing model quality.
Cora, Citeseer, and Pubmed~\cite{sen2008} are citation networks; nodes represent papers connected by citation edges; features are bag-of-word abstracts, and labels represent paper topics.
Amazon PC\xspace and Amazon Photo\xspace~\cite{shchur2018} are subsets of the Amazon co-purchase graph for the computers and photo sections of the website, where nodes represent goods with edges between ones frequently purchased together; node features are bag-of-word reviews, and class labels are product category.
Coauthor CS\xspace and Coauthor PHY\xspace~\cite{shchur2018} are co-authorship networks based on the Microsoft Academic Graph for the computer science and physics fields respectively; nodes are authors, which are connected by edge if they co-authored a paper together; node features are a collection of paper keywords for author’s papers; class
labels indicate most common fields of study.

\spara{Baselines.}
As we study how to leverage the information from both the graph and attributes, we employ two baselines that employ either strictly graph or attribute information.
We give a brief description of the baselines used below:
\vspace{-2mm}
\begin{itemize}[leftmargin=0.2cm,itemindent=.3cm,labelwidth=\itemindent,labelsep=0cm,align=left,topsep=0pt,itemsep=-1ex,partopsep=0ex,parsep=1ex]
    \item \textbf{k-means(features)} is our baseline that only considers the feature data. We use the local Lloyd algorithm~\cite{lloyd1982} with the k-means++ seeding~\cite{arthur2007} strategy.
    \item \textbf{SBM}~\cite{peixoto2014} is a baseline that only relies on the graph structure. We estimate a constrained stochastic block model with given number of $k$ clusters, maximizing the modularity~\cite{newman2006} of the network.
    \item \textbf{k-means(DGI)}~\cite{velickovic2018} demonstrates the need of joint learning of clusters and representations. We learn unsupervised node representation with DGI and run k-means on the resulting representations.
    \item \textbf{MinCutPool(graph, features)}~\cite{bianchi2020} is a deep pooling method that we re-interpret as clustering.
\end{itemize}

\begin{table}[!t]
\small
\centering{
\newcolumntype{R}{>{\raggedleft\arraybackslash}X}
\newcolumntype{C}{>{\centering\arraybackslash}X}
\caption{Results on four datasets from~\cite{sen2008} in terms of graph conductance $\mathcal{C}$, modularity $\mathcal{Q}$, NMI with ground-truth labels, and pairwise F1 measure. We highlight best neural method performance.}
\begin{tabularx}{\linewidth}{p{2cm}CCCCCCCCCCCC@{}}
\toprule
\multicolumn{1}{C}{} & \multicolumn{4}{c}{\textbf{Cora\xspace}} & \multicolumn{4}{c}{\textbf{Citeseer\xspace}} & \multicolumn{4}{c}{\textbf{Pubmed\xspace}} \\
\multicolumn{1}{C}{} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}}  \\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}\cmidrule(lr){10-13}
\emph{method} & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$ & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$ & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$\\
\midrule
\mbox{k-means(features)} & 61.7 & 19.8 & 18.5 & 27.0 & 60.5 & 30.3 & 24.5 & 29.2 & 55.8 & 33.4 & 19.4 & 24.4 \\
\mbox{SBM} & 15.4 & 77.3 & 36.2 & 30.2 & 14.2 & 78.1 & 15.3 & 19.1 & 39.0 & 53.5 & 16.4 & 16.7 \\
\midrule
\mbox{k-means(DGI)} & 28.0 & 64.0 & \textbf{52.7} & 40.1 & 17.5 & 73.7 & \textbf{40.4} & 39.4 & 82.9 & 9.6 & 22.0 & 26.4 \\
\mbox{MinCut} & 23.3 & 70.3 & 35.8 & 25.0 & 14.1 & 78.9 & 25.9 & 20.1 & 29.6 & 63.1 & 25.4 & 15.8 \\
\midrule
\mbox{\textsc{DMoN}\xspace} & \textbf{12.2} & \textbf{76.5} & 48.8 & \textbf{48.8} & \textbf{5.1} & \textbf{79.3} & 33.7 & \textbf{43.2} & \textbf{17.7} & \textbf{65.4} & \textbf{29.8} & \textbf{33.9} \\
\bottomrule
\end{tabularx}}
\label{tbl:clustering}
\end{table}
\begin{table}[!t]
\small
\centering{
\newcolumntype{R}{>{\raggedleft\arraybackslash}X}
\newcolumntype{C}{>{\centering\arraybackslash}X}
\caption{Results on four datasets from~\cite{shchur2018} in terms of graph conductance $\mathcal{C}$, modularity $\mathcal{Q}$, NMI with ground-truth labels, and pairwise F1 measure. We highlight best neural method performance.}
\begin{tabularx}{\linewidth}{@{}p{1cm}CCCCCCCCCCCCCCCC@{}}
\toprule
\multicolumn{1}{C}{} & \multicolumn{4}{c}{\textbf{Amazon PC\xspace}} & \multicolumn{4}{c}{\textbf{Amazon Photo\xspace}} & \multicolumn{4}{c}{\textbf{Coauthor CS\xspace}} & \multicolumn{4}{c}{\textbf{Coauthor PHY\xspace}} \\
\multicolumn{1}{C}{} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}} & \multicolumn{2}{c}{\emph{graph}} & \multicolumn{2}{c}{\emph{labels}}  \\
\cmidrule(lr){2-5}\cmidrule(lr){6-9}\cmidrule(lr){10-13}\cmidrule(lr){14-17}
\emph{method} & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$ & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$ & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$ & $\mathcal{C}\downarrow$ & $\mathcal{Q}\uparrow$ & NMI$\uparrow$ & F1$\uparrow$\\
\midrule
\mbox{k-m(feat)} & 84.5 & 5.4 & 21.1 & 19.2 & 79.6 & 10.5 & 28.8 & 19.5 & 49.1 & 23.1 & 35.7 & 39.4 & 57.0 & 19.4 & 30.6 & 42.9 \\
\mbox{SBM} & 31.0 & 60.8 & 48.4 & 34.6 & 18.6 & 72.7 & 59.3 & 47.4 & 20.3 & 72.7 & 58.0 & 47.7 & 25.9 & 66.9 & 45.4 & 30.4 \\
\midrule
\mbox{k-m(DGI)} & 61.9 & 22.8 & 22.6 & 15.0 & 51.5 & 35.1 & 33.4 & 23.6 & 35.1 & 57.8 & 64.6 & 51.9 & 38.6 & 51.2 & 51.0 & 30.6 \\
\mbox{MinCut} & --- & --- & --- & --- & --- & --- & --- & --- & 22.7 & 70.5 & 64.6 & 47.8 & 27.8 & 64.3 & 48.3 & 24.9 \\
\midrule
\mbox{\textsc{DMoN}\xspace} & \textbf{18.0} & \textbf{59.0} & \textbf{49.3} & \textbf{45.4} & \textbf{12.7} & \textbf{70.1} & \textbf{63.3} & \textbf{61.0} & \textbf{17.5} & \textbf{72.4} & \textbf{69.1} & \textbf{59.8} & \textbf{18.8} & \textbf{65.8} & \textbf{51.9} & \textbf{37.0} \\
\bottomrule
\end{tabularx}}
\label{tbl:clustering2}
\end{table}
\spara{Metrics.}
We measure both the graph-based metrics of clustering and label correlation to study clustering performance of attributed graphs both in terms of graph and attribute structure.
For graph-level metrics, we report average cluster conductance (as per definition from~\cite{yang2015}) and graph modularity~\cite{newman2006}.
For ground-truth label correlation analysis, we report normalized mutual information (NMI) between the cluster assignments and labels and pairwise F-1 score between all node pairs and their associated cluster pairs.
Where possible, we normalize all metrics by multiplying them by 100 for ease of comparison.

\spara{Parameter settings.}
We run all experiments for 10 times and average results across runs.
All models were implemented in Tensorflow 2 and trained on CPUs.
We fix the architecture for all clustering networks to have one hidden layer with 512 neurons in it for real-world datasets and 64 for small synthetic graphs.
We set the number of clusters to 16 for all datasets and methods.

\begin{table}[!b]
\vspace{-2em}
\small
\centering
\newcolumntype{S}{>{\hsize=.1\hsize}X}
\newcolumntype{L}{>{\hsize=.3\hsize}X}
\caption{\label{tab:synthetic-scenarios} Synthetic ADC-SBM benchmark scenarios.}
\begin{tabularx}{\textwidth}{SLX}
\toprule
Scenario & Parameter & Description \\
\midrule
1 & $d_{out}\in[2.0, 5.0]$ & Increase graph cluster mixing signal. Higher = weaker clusters. \\
2 & $\sigma_c\in[10^{-2}, 10^1]$ & Increase feature cluster center variance. Higher = stronger clusters.\\
3 & $\sigma_c\in[10^{-2}, 10^1]$ & Increase feature cluster center variance, with nested feature clusters.\\
4 & $\sigma_c\in[10^{-2}, 10^1]$ & Increase feature cluster center variance, with grouped feature clusters.\\
5 & $d\in [2^2, 2^7]$ & Increase average degree. Higher = clearer graph signal.\\
6 & $d_{max} \in [2^2, 2^{10}]$ & Increase power law upper-bound. Higher = more extreme power law.\\
\bottomrule
    \end{tabularx}
\end{table}
\subsection{Simulation Experiments on Stochastic Block Model}\label{ssec:sbm-experiments}
To better explore the robustness of our approach to variance in the graph and node features, we propose a study on synthetic graphs using an \emph{attributed}, \emph{degree-corrected} stochastic block model (ADC-SBM).
The SBM~\cite{snijders1997estimation} plants a partition of clusters (``blocks'') in a graph, and generates edges via a distribution conditional on that partition.
This model has been used extensively to benchmark graph clustering methods~\cite{fortunato2016community}, and has recently been used for experiments on state-of-the-art GNNs~\cite{dwivedi2020benchmarking}.
In our version of the model, node features are also generated, using a multivariate mixture model, with the mixture memberships having some correspondence to the cluster memberships.

To generate an instance of the ADC-SBM, we fix a number of nodes $n$ and a number of clusters $k$, and choose node cluster memberships uniformly-at-random. Define the matrix $\mathbf{D}_{k\times k}$ where $\mathbf{D}_{ij}$ is the expected number of edges shared between nodes in clusters $i$ and $j$. We determine $\mathbf{D}$ by fixing (1) the expected \emph{average} degree of the nodes $d\in \{1, n\}$, and (2) the expected \emph{average} sub-degree $d_{out} \leq d$ of a node to any cluster other than its own.
Note that the difference $d_{in} - d_{out}$, where $d_{in}:=d - d_{out}$, controls the spectral detectability of the clusters~\cite{nadakuditi2012graph}.
Finally, we generate a power-law $n$-vector $\theta$, where $\theta_i$ is proportional to $i$'s expected degree. With the memberships and the parameters $D$ and $\theta$. We use \texttt{graph-tool}~\cite{peixoto_graph-tool_2014} generate the graphs.

\begin{figure}[!t]
\begin{subfigure}[b]{0.24\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/adj.png}%
\caption{}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.24\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/cov1.png}%
\caption{}
\end{subfigure}\hfill
\begin{subfigure}[b]{.24\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/cov2.png}%
\caption{}
\end{subfigure}\hfill
\begin{subfigure}[b]{0.24\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/cov3.png}
\caption{}
\end{subfigure}
\caption{\label{fig:synthetic-illustration}\textbf{Illustration of synthetic data.} a): 4-cluster graph adjacency matrix. b): Covariance matrix of ``matched'' features: features that are clustered according to the graph clusters. c): Covariance matrix of ``nested'' features: features that are clustered by nesting of the graph clusters. d): Covariance matrix of ``grouped'' features: features that are clustered by grouping of the graph clusters.}
\vspace{-1.5em}
\end{figure}
To generate $s$-dimensional features, we first generate feature memberships from $k_f$ cluster labels.
For graph clustering GNNs that operate both on edges and node features, it is important to examine performance on data where feature clusters diverge from or segment the graph clusters: thus potentially $k_f\ne k$. We examine cases where feature memberships \emph{match}, \emph{group}, or \emph{nest} the graph memberships, as illustrated in Figure~\ref{fig:synthetic-illustration}.
With feature memberships in-hand, we generate $k$ zero-mean feature cluster centers from a $s$-multivariate normal with covariance matrix $\sigma_c^2\cdot \mathbf{I}_{s\times s}$.
Then, for feature cluster $i\leq k_f$, we generate its features from a $s$-multivariate normal with covariance matrix $\sigma^2\cdot \mathbf{I}_{s\times s}$.
Note that the ratio $\sigma^2_c/\sigma^2$ controls the expected value of the classical between/within-sum-of-squares of the clusters.

The above paragraphs describe a single generation of our synthetic benchmark model, the ADC-SBM. To explore robustness, we define a ``default'' ADC-SBM, and explore model parameters in a range around the defaults.
We configure our default model as follows: we generate graphs with $n=1,000$ nodes grouped in $k=4$ clusters, and $s=32$-dimensional features grouped in $k_f=4$ matching feature clusters with $\sigma=1$ intra-cluster center variance and $\sigma_c=3$ cluster center variance. We try to model real-world graphs' degree distribution with $d=20$ average degree and $d_{out}=2$ average inter-cluster degree with power law parameters $d_{min}=2, d_{max}=4, \alpha=2$.
In total, we consider 6 different scenarios, as described in Table~\ref{tab:synthetic-scenarios}.

\begin{wrapfigure}[10]{R}{0.3\linewidth}
\vspace{-6mm}
\centering
\begin{subfigure}[b]{0.475\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/mincut.png}%
\caption{}
\end{subfigure}\hfill
\begin{subfigure}[b]{.475\linewidth}
\centering
\includegraphics[width=\linewidth]{figures/mgcn2.png}%
\caption{}
\end{subfigure}
\vspace*{-3mm}
\caption{\label{fig:viz}\small(a) MinCutPool fails to recover ground-truth clusters while (b) \textsc{DMoN}\xspace perfectly recovers ground-truth partition.}
\end{wrapfigure}
Figure~\ref{fig:synthetic-results} demonstrates overwhelming superiority of \textsc{DMoN}\xspace{} in all considered scenarios.
In scenario 1 we see that \textsc{DMoN}\xspace{} can effectively leverage the feature signal to obtain outstanding clustering performance even when the graph structure is close to random, far beyond the spectral detectability threshold (pictured in gray).
Scenario 2 demonstrates that even in the presence of a weak feature signal \textsc{DMoN}\xspace outperforms stochastic SBM minimization, while MinCutPool needs \emph{two orders of magnitude} stronger signal to attain the same performance.
Scenario 3 demonstrates that \textsc{DMoN}\xspace{} is not susceptible to nested feature clusters and can correctly recover underlying graph structure.
Scenario 4 shows that \textsc{DMoN}\xspace is more susceptible to grouped features, as it becomes hard to differentiate graph clusters with similar features, however, MinCutPool is struggling to pick up \emph{any} signal in the data.
Scenarios 5\&6 demonstrate that \textsc{DMoN}\xspace is stable in terms of the degree distribution of the graph, while MinCutPool is not.
We also notice that while the k-means(DGI) baseline offers some improvements over using features or the graph structure alone, it never surpasses the strongest signal provider in the graph, never being better than the best one between k-means(features) and SBM.
Finally, we provide an illustration of the how MinCutPool collapsing into a single cluster while \textsc{DMoN}\xspace succeeds in Figure~\ref{fig:viz}.

\begin{figure}[!t]
\centering
\begin{tikzpicture}
\begin{groupplot}[group style={
                      group name=myplot,
                      group size= 3 by 2, horizontal sep=0.75cm,vertical sep=1.3cm},height=3.75cm,width=0.375 \linewidth,ymin=0,ymax=100,title style={at={(0.5,0.9)},anchor=south},every axis x label/.style={at={(axis description cs:0.5,-0.15)},anchor=north}]
\nextgroupplot[
 	title = \textbf{Scenario 1},
 	legend columns=5,
	legend style={at={(1.675,1.25)},anchor=south},
    legend entries={\textsc{DMoN}\xspace{}, SBM, k-means(features), MinCut, k-means(DGI)},
	ylabel=NMI$\times100$,
	xlabel=$d_{out}$,
	xmin=1,
	xmax=5
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario1.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario1.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario1.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario1.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario1.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario1.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario1.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario1.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario1.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario1.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario1.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario1.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario1.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario1.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario1.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];

\draw [color=cyclegray, ultra thick, draw opacity=0.75] (2.76393202250021,0) -- (2.76393202250021,100);
\draw[->,thick,color=black](axis cs:3.15,28.5)--(axis cs:2.85,60);
\node[anchor=east] (text) at (axis cs:3.9,20){\tiny \textbf{detectability limit}};
\nextgroupplot[
 	title = \textbf{Scenario 2},
	yticklabels={,,},
	xlabel=$\sigma_c$,
	xmode=log,
	xmin=0.01,
	xmax=10
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario2.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario2.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario2.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario2.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario2.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario2.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario2.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario2.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario2.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario2.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario2.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario2.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario2.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario2.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario2.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];
\nextgroupplot[
 	title = \textbf{Scenario 3},
	yticklabels={,,},
	xlabel=$\sigma_c$,
	xmode=log,
	xmin=0.01,
	xmax=10
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario3.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario3.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario3.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario3.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario3.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario3.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario3.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario3.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario3.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario3.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario3.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario3.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario3.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario3.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario3.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];
\nextgroupplot[
 	title = \textbf{Scenario 4},
	xlabel=$\sigma_c$,
	ylabel=NMI$\times100$,
	xmode=log,
	xmin=0.01,
	xmax=10
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario4.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario4.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario4.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario4.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario4.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario4.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario4.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario4.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario4.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario4.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario4.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario4.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario4.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario4.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario4.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];
\nextgroupplot[
 	title = \textbf{Scenario 5},
	yticklabels={,,},
	xlabel=$d$,
	xmode=log,
	xmin=4,
	xmax=128,
	log basis x={2}
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario5.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario5.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario5.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario5.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario5.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario5.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario5.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario5.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario5.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario5.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario5.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario5.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario5.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario5.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario5.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];
\nextgroupplot[
 	title = \textbf{Scenario 6},
	yticklabels={,,},
	xlabel=$d_{\max}$,
	xmode=log,
	xmin=4,
	xmax=1024,
	log basis x={2},
    max space between ticks=20,
]
\addplot[very thick,color=cycle2] table[x=param,y=mgcn2] {data/scenario6.tex};
\addplot[very thick,color=cycle3] table[x=param,y=sbm] {data/scenario6.tex};
\addplot[very thick,color=cycle4] table[x=param,y=km] {data/scenario6.tex};
\addplot[very thick,color=cycle5] table[x=param,y=mincut] {data/scenario6.tex};
\addplot[very thick,color=cycle6] table[x=param,y=dgi] {data/scenario6.tex};

\addplot[name path=mgcn_top_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_high] {data/scenario6.tex};
\addplot[name path=mgcn_btm_,color=cycle2!70,forget plot] table[x=param,y=mgcn2_low] {data/scenario6.tex};
\addplot[cycle2!50,fill opacity=0.5,forget plot] fill between[of=mgcn_top_ and mgcn_btm_];

\addplot[name path=sbm_top_,color=cycle3!70,forget plot] table[x=param,y=sbm_high] {data/scenario6.tex};
\addplot[name path=sbm_btm_,color=cycle3!70,forget plot] table[x=param,y=sbm_low] {data/scenario6.tex};
\addplot[cycle3!50,fill opacity=0.5,forget plot] fill between[of=sbm_top_ and sbm_btm_];

\addplot[name path=km_top_,color=cycle4!70,forget plot] table[x=param,y=km_high] {data/scenario6.tex};
\addplot[name path=km_btm_,color=cycle4!70,forget plot] table[x=param,y=km_low] {data/scenario6.tex};
\addplot[cycle4!50,fill opacity=0.5,forget plot] fill between[of=km_top_ and km_btm_];

\addplot[name path=mincut_top_,color=cycle5!70,forget plot] table[x=param,y=mincut_high] {data/scenario6.tex};
\addplot[name path=mincut_btm_,color=cycle5!70,forget plot] table[x=param,y=mincut_low] {data/scenario6.tex};
\addplot[cycle5!50,fill opacity=0.5,forget plot] fill between[of=mincut_top_ and mincut_btm_];

\addplot[name path=dgi_top_,color=cycle6!70,forget plot] table[x=param,y=dgi_high] {data/scenario6.tex};
\addplot[name path=dgi_btm_,color=cycle6!70,forget plot] table[x=param,y=dgi_low] {data/scenario6.tex};
\addplot[cycle6!50,fill opacity=0.5,forget plot] fill between[of=dgi_top_ and dgi_btm_];
\end{groupplot}
\end{tikzpicture}
\vspace*{-5mm}
\caption{Synthetic results on the ADC-SBM model with 6 different scenarios described in Table~\ref{tab:synthetic-scenarios}.}\label{fig:synthetic-results}
\end{figure}

\subsection{Real-world Data}\label{ssec:real-experiments}
We now move on to studies on real-world networks, featuring \textsc{DMoN}\xspace and 4 baselines on 7 different datasets.
\textsc{DMoN}\xspace achieves better clustering performance than its neural counterparts on every single dataset and metric besides losing twice to DGI+k-means on Cora and Citeseer in terms of NMI.
Compared to SBM, a method that exclusively optimizes modularity, we are able to stay within 2 NMI percentage points, while simultaneously clustering the features.
Surprisingly, on Citeseer and Pubmed we achieve better modularity than the method designed to optimize it -- we attribute that to very high correlation between graph structure and the features.
We also highlight that we beat MinCutPool in terms of conductance (average graph cut) on all datasets, even though it attempts to optimize for this metric.

On Amazon PC\xspace and Amazon Photo\xspace MinCutPool failed to converge, even with tuning the orthogonality regularization parameter.
We attribute that to extremely uneven structure of these graphs, as popular products are co-purchased with a lot of other items, so the effects discussed in~\cite{epasto2017ego,leskovec2008statistical} are prohibiting good cuts.
This corresponds to high values of $d$ ad $d_{\max}$ in our synthetic scenarios 5 and 6.

Overall, \textsc{DMoN}\xspace demonstrates excellent performance on both graph clustering and label correlation, successfully leveraging both graph and attribute information.
Both synthetic and real-world experiments prove that \textsc{DMoN}\xspace is vastly superior to its counterparts in attributed graph clustering.

 %
\section{Conclusion}\label{sec:conclusion}
In this work, we study GNN pooling through the lens of attributed graph clustering.
We introduce Deep Modularity Networks (\textsc{DMoN}\xspace), an unsupervised objective and realize it with a GNN which can recover high quality clusters.
We compare against challenging baselines that baselines that optimize structure (SBM), features (kmeans), or both (DGI+k-means), in addition to a recently proposed state-of-the-art pooling method (MinCutPool).

We explore the limits of GNN clustering methods in terms of both graph and feature signals using synthetic data, where we see that \textsc{DMoN}\xspace better leverages structure and attributes than all existing methods.
In extensive experiments on real datasets we show that the clusters found by \textsc{DMoN}\xspace are more likely to correspond to ground truth labels, and have better properties as illustrated by clustering metrics (e.g.\ conductance or modularity).
We hope that this work will further advancements in unsupervised learning for GNNs as well as attributed graph clustering, allowing further advances in graph learning.

\newpage

\section*{Broader Impact}
We believe that better understanding of graph pooling through the lens of attributed graph clustering will allow for the design of more effective graph neural network architectures. Our reformulation encourages explicit grouping of nodes in graphs, which we envision will be useful in domains with naturally occurring community structure, e.g.\ social networks or information retrieval. Additionally, our M-GCN could serve as the ``graph-context'' head in semi-supervised architectures \cite{yang2016revisiting}, allowing for stronger community-based constraints in those systems.

Our work also allows translation of open questions in network science about attributed graph clustering to the GNN domain, potentially enabling powerful new modelling solutions. Specifically, existing combinatorial \cite{yang2013community} and bayesian \cite{newman2016structure} methods for attributed graph clustering use graph features to improve clustering accuracy. \cite{newman2016structure} shows that using features allows their estimation-based approach to overcome the spectral ``detectability'' limit for graphs \cite{nadakuditi2012graph}, a result that we replicate in this study. Using our model, more ideas from the network science community may be ported to GNN architectures. \bibliographystyle{plain}