# Modularity aided consistent attributed graph clustering via coarsening

Samarth Bhatia*<EMAIL_ADDRESS>
Indian Institute of Technology, Delhi Yukti Makhija*
<EMAIL_ADDRESS>
Indian Institute of Technology, Delhi Manoj Kumar<EMAIL_ADDRESS>
Indian Institute of Technology, Delhi Sandeep Kumar<EMAIL_ADDRESS>
Indian Institute of Technology, Delhi

* denotes equal contribution

###### Abstract

Graph clustering is an important unsupervised learning technique for
partitioning graphs with attributes and detecting communities. However,
current methods struggle to accurately capture true community structures and
intra-cluster relations, be computationally efficient, and identify smaller
communities. We address these challenges by integrating coarsening and
modularity maximization, effectively leveraging both adjacency and node
features to enhance clustering accuracy. We propose a loss function
incorporating log-determinant, smoothness, and modularity components using a
block majorization-minimization technique, resulting in superior clustering
outcomes. The method is theoretically consistent under the Degree-Corrected
Stochastic Block Model (DC-SBM), ensuring asymptotic error-free performance
and complete label recovery. Our provably convergent and time-efficient
algorithm seamlessly integrates with graph neural networks (GNNs) and
variational graph autoencoders (VGAEs) to learn enhanced node features and
deliver exceptional clustering performance. Extensive experiments on benchmark
datasets demonstrate its superiority over existing state-of-the-art methods
for both attributed and non-attributed graphs.

##  1 Introduction

Clustering is an unsupervised learning method that groups nodes together based
on their attributes or graph structure, without considering node labels. This
versatile approach has numerous applications in diverse fields, such as social
network analysis (Tang et al., 2008), genetics, and bio-medicine (Cheng & Ma,
2022; Buterez et al., 2021), knowledge graphs (Hamaguchi et al., 2017), and
computer vision (Mondal et al., 2021; Caron et al., 2018). The wide range of
applications has led to the development of numerous graph clustering
algorithms designed for specific challenges within these domains. State-of-
the-art graph clustering methods predominantly fall into cut-based,
similarity-driven, or modularity-based categories.

Cut-based methods (Wei & Cheng, 1989; Shi & Malik, 2000; Bianchi et al.,
2020), aiming to minimize the number of edges (or similar metric) in a cut,
may fall short in capturing the true community structure if the cut’s edge
count doesn’t significantly deviate from random graph expectations. This
approach originated from the Fiedler vector, which yields a graph cut with a
minimal number of edges over all possible cuts (Fiedler, 1973; Newman, 2006a).
Similarity-based techniques, reliant on pairwise node similarities for
clustering leverage metrics of similarity to group nodes with shared
characteristics. These can be computationally intensive and susceptible to
noise, potentially yielding suboptimal results, especially in sparse or noisy
data scenarios. Modularity-based methods rely on a statistical approach,
measuring the disparity in edge density between a graph and a random graph
with the same degree sequence. On the other hand, these methods exhibit a
resolution limit (Fortunato & Barthélemy, 2007), as the emphasis on maximizing
modularity may inadvertently lead to the neglect of smaller community
structures within the graph. Each clustering approach, whether cut-based,
similarity-based, or modularity-based, bears its own set of limitations that
necessitates careful consideration based on the characteristics of the
underlying graph data.

Moreover, graph reduction techniques such as coarsening, summarization, or
condensation can also be utilized to facilitate the clustering task (Dhillon
et al., 2007; Kumar et al., 2023; Loukas & Vandergheynst, 2018; Loukas, 2019).
In the context of graph coarsening, the objective is to learn a reduced graph
by merging similar nodes into supernodes. While this coarsening process is
traditionally employed for graph reduction, it can be extended to clustering
by reducing the original graph such that each class corresponds to a
supernode. However, relying solely on the loss function of coarsening may lead
to suboptimal performance, particularly when reducing the graph size to the
number of classes results in a considerable reduction and may lead to
significant loss of information.

In this study, we introduce an optimization-based framework designed to
enhance the clustering task by leveraging both the adjacency and the node
features of the graph. Our proposed framework strategically incorporates
coarsening and modularity maximization, refining partitioning outcomes and
bolstering the overall effectiveness of the clustering process. The algorithm
minimizes a nuanced loss function, encompassing a $\log\det$ term, smoothness,
and modularity components, to ensure efficient clustering. Formulated as a
multi-block non-convex optimization problem, our approach is adeptly addressed
through a block majorization-minimization technique, wherein variables are
updated individually while keeping the remainder fixed. The resulting
algorithm demonstrates convergence, showcasing its efficacy in efficiently
solving the proposed optimization problem.

To enhance the clustering performance, we embed our Q-MAGC objective function
into various Graph Neural Network (GNN) architectures, introducing the Q-GCN
algorithm. This novel approach elevates the quality of learned representations
by leveraging the message passing and aggregation mechanisms of Graph
Convolutional Networks (GCNs), ultimately improving the clustering outcomes.
An additional feature of our technique is its capacity to explore inter-
cluster relationships. The node attributes derived from the coarsened graph at
the conclusion of the process serve as cluster embeddings, shedding light on
the distinct characteristics of each cluster. Simultaneously, the edges within
the coarsened graph unveil valuable insights into the relationships and
connections between different clusters, providing a comprehensive
understanding of the overall graph structure. This is a contributing factor to
the improvement observed over existing methods and is particularly significant
for conducting first-hand analyses on large unlabelled datasets.

We introduce two additional algorithms, Q-VGAE and Q-GMM-VGAE, incorporating
variational graph auto-encoders to further enhance clustering accuracy.
Through comprehensive experiments, we demonstrate the effectiveness of our
proposed algorithms, clearly surpassing the performance of current state-of-
the-art methods on benchmark datasets. Our innovative approach, showcases
notable improvements in clustering performance, thereby solidifying the
robustness and superiority of our proposed methods in handling real-world
benchmark datasets.

Key Contributions.

  * •

We present the first optimization-based framework for attributed graph
clustering through coarsening via modularity maximization. Our approach
demonstrates efficiency, theoretical convergence, and addresses limitations of
existing methods. The paper offers comprehensive analysis and provides
theoretical guarantees including KKT optimality, and convergence analysis
which are often absent in prior research.

  * •

We show that our method is theoretically (weakly and strongly) consistent
under a Degree-Corrected SBM (DC-SBM) and shows asymptotically no errors
(weakly consistent) and complete recovery of the original labels (strongly
consistent).

  * •

We demonstrate the seamless integration of the proposed clustering objective
with GNN-based architectures, leveraging message-passing in GNNs to enhance
the performance of our method. This is also backed up by experiments. 111Refer
to Appendix K for the code

  * •

We perform thorough experimental validation on a diverse range of real-world
and synthetic datasets, encompassing both attributed and non-attributed graphs
of varying sizes. The results demonstrate the superior performance of our
method compared to existing state-of-the-art approaches. We want to highlight
that our method does not specialise for very large graphs, yet we present
preliminary results in this regard.

  * •

We conduct ablation studies to evaluate the behavior of the loss terms,
compare runtime and complexities, and perform a comprehensive evaluation of
modularity.

Notations. Let $G=\\{V,E,A,X\\}$ be a graph with node set
$V=\\{v_{1},v_{2},...,v_{p}\\}\ (|V|=p)$, edge set $E\subset V\times
V\\}(|E|=e)$, weight (adjacency) matrix $A$ and node feature matrix
$X\in\mathbb{R}^{p\times n}$. Also, let
$\textbf{d}=A.\textbf{1}\in\mathbb{Z}_{+}^{p}$ be the degree vector, where 1
is the vector pf size $p$ having all entries 1. Then, the graph Laplacian is
$\Theta=\textrm{diag}(\textbf{d})-A$ and the set of all valid Laplacian
matrices is defined as: $S_{\Theta}=\\{\Theta\in\mathbb{R}^{p\times
p}|\Theta_{ij}=\Theta{ji}\leq 0\ \text{for}\ i\neq
j,\Theta_{ii}=\sum_{j=1}^{p}\Theta_{ij}\\}$.

##  2 Related Works

In this section, we review relevant existing works and highlight their
limitations, thereby motivating the need for an improved clustering
formulation.

Graph Clustering via Coarsening. Graph coarsening can be extended to graph
clustering by reducing the size of the coarsened graph to the number of
classes, i.e., $k=l$. However, in most graphs, the number of classes is very
small, and reducing the graph to this extent may lead to poor clustering
quality. For example, DiffPool (Ying et al., 2018) learns soft cluster
assignments at each layer of the GNN and two extra losses, entropy to penalize
the soft assignments and a link prediction based loss. Next, SAGPool (Lee et
al., 2019) calculates attention scores and node embeddings to determine the
nodes that need to be preserved or removed. Top-k (Gao & Ji, 2019) also works
by sparsifying the graph with the learned weights. MinCutPool (Bianchi et al.,
2020) formulates a differentiable relaxation of spectral clustering via
pooling. However, (Tsitsulin et al., 2023) show that MinCutPool’s orthogonal
regularization dominates over the clustering objective and the objective is
not optimized. Some disadvantages of these methods are stability and
computational complexity in the case of SAGPool and DiffPool and convergence
in MinCutPool. To address these challenges, we need a better loss function to
perform the clustering task effectively.

Deep Graph Clustering. Previous literature can be classified based on
contrastive and non-contrastive methods. On the non-contrastive side, Pan et
al. (2018) proposed ARGA and ARVGA, enforcing the latent representations to
align to a prior using adversarial learning. By utilizing an attention-based
graph encoder and a clustering alignment loss, Wang et al. (2019) propose
DAEGC. Liu et al. (2022) design the DCRN model to alleviate representation
collapse by a propagation regularization term minimizing the JSD between the
latent and its product with normalized $A$. Contrastive methods include AGE
(Cui et al., 2020) which builds a training set by adaptively selecting node
pairs that are highly similar or dissimilar after filtering out high-frequency
noises using Laplacian smoothing. Zhao et al. (2021) propose GDCL to correct
the sampling bias by choosing negative samples based on the clustering label.

VGAEs (Kipf & Welling, 2016a) are an increasingly popular class of GNNs that
leverage variational inferenceKingma & Welling (2014) for learning latent
graph representations in unsupervised settings. They reconstruct the adjacency
matrix after passing the graph through an encoder-decoder architecture. Many
attempts have been made to use VGAEs with k-means on latent embeddings, but it
has been unsuitable for clustering. This is primarily because embedded
manifolds obtained from VGAEs are curved and must be flattened before any
clustering algorithms using Euclidean distance are applied. Refer to Appendix
Section L for a detailed explanation. VGAEs only use a single Gaussian prior
for the latent space, whereas clustering requires the integration of meta-
priors. Additionally, the inner-product decoder fails to capture locality and
cluster information in the formed edges. Several clustering-oriented variants
of VGAEs (Mrabah et al., 2022; Hui et al., 2020) have been developed that
overcome most of these challenges. GMM-VGAE (Hui et al., 2020) partitions the
latent space using a Gaussian Mixture Model and assigns a separate prior for
each cluster to better model complex data distributions. Despite the
improvement in performance, it’s inner-product decoder cannot capture locality
information.

Modularity Maximization. Various heuristic algorithms have been established
that solve the NP-hard problem of modularity maximization including sampling,
simulated annealing (Guimerà & Nunes Amaral, 2005; Guimerà & Amaral, 2005),
and greedy algorithms (Louvain/Leiden) (Newman, 2004; Blondel et al., 2008a).
These algorithms require intensive compute and don’t use node features.
Modularity maximization using GNNs has also garnered attention recently. DMoN
(Tsitsulin et al., 2023) optimizes only for modularity (with a collapse
regularization to prevent the trivial solution) but offers no theoretical
guarantees about convergence. Modularity-Aware GAEs and VGAEs (Salha-Galvan et
al., 2022) use a prior membership matrix using Louvain algorithm and optimize
for modularity using RBF kernel as a same-community assignment proxy.
DGCLUSTER (Bhowmick et al., 2023) is a semi-supervised method that makes use
of either a subset of labels or pairwise memberships as part of their
"Auxiliary Information" and a modularity objective.

##  3 Background

In this section, we introduce the concepts that play a key role in the
formulation of our method.

###  3.1 Graph Coarsening

Graph coarsening is a graph dimensionality reduction technique used in large-
scale machine learning to construct a smaller or coarsened graph $G_{c}$ from
the original graph $G=\\{V,E,A,X\\}$ while preserving properties of the
original graph $G$. Given an original graph $G(V,E,A,X)$, graph coarsening
aims to learn a mapping matrix $C\in\mathbb{R}_{+}^{p\times k}$, where $p$ is
the number of nodes in the original graph and $k$ is the number of nodes in
the coarsened graph. Each non-zero entry of the mapping matrix $C$, i.e.,
$C_{ij}$, indicates that the $i$-th node of the original graph is mapped to
the $j$-th supernode. Moreover, for a balanced mapping, the mapping matrix $C$
belongs to the following set Kumar et al. (2023); Loukas & Vandergheynst
(2018); Loukas (2019):

$\displaystyle\mathcal{C}=\Big{\\{}C\in\mathbb{R}_{+}^{p\times k}|\ \langle
C_{i},C_{j}\rangle=0\ $ $\displaystyle\forall\;i\neq j,\quad\langle
C_{i},C_{i}\rangle=d_{i},\left\|C_{i}\right\|_{0}\geq 1\ \text{and}\
\left\|[C^{\top}]_{i}\right\|_{0}=1\Big{\\}}$ (1)

For $C\in\mathcal{C}$, the relationship between the original graph Laplacian
$\Theta$, the coarsened graph Laplacian matrix $\Theta_{c}$, and the mapping
matrix $C$ is given by $\Theta_{c}=C^{T}\Theta C$.

###  3.2 Spectral Modularity Maximization

Spectral Clustering is the most direct approach to graph clustering, where we
minimize the volume of inter-cluster edges (i.e., the total number of edges
between clusters, not within them). Modularity, introduced in (Newman, 2006b),
is the difference between the number of edges within a cluster
$\mathbf{C_{i}}$ and the expected number of such edges in a random graph with
an identical degree sequence. It is mathematically defined as:

$\displaystyle\mathcal{Q}$
$\displaystyle=\frac{1}{2e}\sum_{i,j=1}^{k}\Big{[}A_{ij}-\frac{d_{i}d_{j}}{2e}\Big{]}\delta(c_{i},c_{j})$
(2)

where $\delta(c_{i},c_{j})$ is the Kronecker delta between clusters $i$ and
$j$. Maximizing this form of modularity is NP-hard (Brandes et al., 2008).
However, we can approximate it using a spectral relaxation, which involves a
modularity matrix $B$. The modularity matrix $B$ and spectral modularity are
defined as follows:

$\displaystyle B$
$\displaystyle=A-\frac{\textbf{dd}^{T}}{2e},\qquad\mathcal{Q}=\frac{1}{2e}Tr(C^{T}BC)$
(3)

Similar to the graph Laplacian matrix $\mathbf{\Theta}$, $\mathbf{B}$ is
symmetric and is defined such that its row-sums and column-sums are zero,
thereby making $\mathbf{1}$ one of its eigenvectors and $0$ the corresponding
eigenvalue. These spectral properties of the modularity matrix are also
observed in the Laplacian, as noted in Newman (2006b), which is a crucial
element in spectral clustering. Modularity is maximized when $u_{1}^{T}s$ is
maximized, where $u$ are the eigenvectors of $B$ and $s$ is the community
assignment vector, i.e., placing the majority of the summation in $Q$ on the
first (and largest) eigenvalue of $B$. Moreover, modularity is closely
associated with community detection. These special spectral properties make
$\mathbf{B}$ an ideal choice for graph clustering. While modularity
maximization has been extensively studied, heuristic algorithms for it are
computationally intensive (such as the Newman-Girvan algorithm with
$\mathcal{O}(n^{3})$ time complexity). The Louvain/Leiden algorithms improve
on this.

###  3.3 Motivation and Proposed Problem Formulation

Current graph clustering methods often struggle to accurately capture both
intra-cluster and inter-cluster relationships, achieve computational
efficiency, and identify smaller communities, thereby limiting their
effectiveness. These methods face significant challenges including
instability, high computational complexity, and the lack of convergence
guarantees, especially when employing coarsening techniques. Modularity
maximization, despite extensive study, relies on computationally intensive
heuristics and lacks theoretical convergence guarantees. To address these
issues, we propose an optimization-driven framework that strategically
integrates coarsening and modularity. By incorporating both adjacency and node
features, our approach aims to robustly capture intra-cluster and inter-
cluster dynamics. This framework ensures stability, guarantees convergence,
and consistently delivers superior clustering results with enhanced
computational efficiency compared to existing methods. Given original graph
$G(V,E,A,X)$, the proposed optimization formulation is:

$\displaystyle\underset{C}{\min}\;\mathcal{L}_{MAGC}$ $\displaystyle=\
f(C,X,\Theta)+g(C,A)+h(C,\Theta)$ $\displaystyle\text{subject to}\ C$
$\displaystyle\in\mathcal{S_{C}}=\\{C\in\mathbb{R}^{p\times l}|C\geq 0,\
\left\|[C^{T}]_{i}\right\|_{2}^{2}\leq 1\\}\;\ \forall\ i=1,2\ldots p$ (4)

Here, $C\in\mathbb{R}_{+}^{p\times l}$ represents the clustering matrix to be
learned, where each non-zero entry $C_{ij}$ indicates that the $i$-th node is
assigned to the $j$-th cluster. The term $f(C,X,\Theta)$ denotes the graph
coarsening objective, managing inter-cluster edges. The function $g(C,A)$
represents the modularity objective, improving clustering performance.
Additionally, $h(C,\Theta)$ is a regularization term enforcing desirable
properties in the clustering matrix, as defined in Equation 1. The overarching
goal of the optimization problem 3.3 is to choose $f(C,X,\Theta)$, $g(C,A)$,
and $h(C,\Theta)$ such that nodes are optimally clustered and inter-cluster
connectivity is maintained. In the next section, we will develop the
clustering algorithm, leveraging the coarsening objective and modularity to
enhance clustering performance.

##  4 Proposed Method

Given a graph $G=\\{V,E,A,X\\}$, considering
$f(C,X_{C},\Theta)=\text{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$,
$g(C,A)=-\frac{\beta}{2e}\textrm{tr}(C^{T}BC)$, and
$h(C,\Theta)=-\gamma\log\det(C^{T}\Theta C+J)$, to obtain clustering matrix
$C$ we formulate the following optimization problem :

$\displaystyle\underset{X_{C},C}{\min}\;\mathcal{L}_{MAGC}$ $\displaystyle=\
\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})-\frac{\beta}{2e}\textrm{tr}(C^{T}BC)-\gamma\log\det(C^{T}\Theta C+J)$
$\displaystyle\text{subject to}\ C\in\mathcal{S_{C}}$
$\displaystyle=\\{C\in\mathbb{R}^{p\times
l}|\left\|C^{T}_{i}\right\|_{2}^{2}\leq 1\\}\;\forall
i,X=CX_{C}\;\;\text{where, }J=\frac{1}{k}\textbf{1}_{k\times k}$ (5)

The term $\textrm{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$ represents the smoothness
or Dirichlet energy of the coarsened graph while $C^{T}\Theta C$ is the
Laplacian matrix of the coarsened graph. Minimzing smoothness ensures that the
clusters or supernodes with similar features are linked with stronger weights.
The term $\textrm{tr}(C^{T}BC)$ corresponds to the graph’s modularity,
enhancing the quality of the clusters formed. The term $-\log\det(C^{T}\Theta
C+J)$ is crucial for maintaining inter-cluster edges in the coarsened graph.
For a connected graph matrix with $l$ super-nodes or clusters, the rank of
$C^{T}\Theta C$ is $l-1$. Adding $J$ to $C^{T}\Theta C$ makes $C^{T}\Theta
C+J$ a full-rank matrix without altering the row and column space of
$C^{T}\Theta C$ (Kumar et al., 2020).

Problem 4 is a multi-block non-convex optimization problem solved using the
Block Successive Upper-bound Minimization (BSUM) framework. All terms except
modularity are convex in nature, which we prove in Appendix A. We iteratively
update $C$ and $X_{C}$ alternately while keeping the other constant. This
process continues until convergence or the stopping criteria are met. Since
the constraint $X=CX_{C}$ is hard and difficult to enforce, we relax it by
adding the term $\frac{\alpha}{2}\|X-CX_{C}\|_{F}^{2}$. This term ensures that
each node is assigned to a cluster, leaving no node unassigned. Additionally,
when needed, we can add an optional sparsity regularization term
$\lambda\|C^{T}\|_{1,2}^{2}$, which can be seen as ensuring that each node is
assigned to exactly one cluster, avoiding any overlap in node assignments
across clusters.

Update rule of C. Treating $X_{C}$ as constant and $C$ as a variable the sub-
problem for $C$ is:

$\displaystyle\underset{C}{\min}\;f(C)$ $\displaystyle=\
\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})-\frac{\beta}{2e}\textrm{tr}(C^{T}BC)-\gamma\log\det(C^{T}\Theta
C+J)+\frac{\alpha}{2}\|X-CX_{C}\|_{F}^{2}$ $\displaystyle\text{subject to}\ C$
$\displaystyle\in\mathcal{S_{C}}=\\{C\in\mathbb{R}^{p\times l}|C\geq 0,\
\left\|C^{T}_{i}\right\|_{2}^{2}\leq 1\\}\;\forall i,\;\;\text{where,
}J=\frac{1}{k}\textbf{1}_{k\times k}$ (6)

By using the first-order Taylor series approximation, a majorised function for
$f(C)$ at $C^{(t)}$ can be obtained as:

$\displaystyle g(C|C^{t})$ $\displaystyle=f(C^{t})+(C-C^{t})\nabla
f(C^{t})+\frac{L}{2}\left\|C-C^{t}\right\|^{2}$ (7)

where $f(C)$ is $L-$Lipschitz continuous gradient function
$L=\max(L_{1},L_{2},L_{3},L_{4})$ with $L_{1},L_{2},L_{3},L_{4}$ the Lipschitz
constants of $-\gamma\text{log det}(C^{T}\Theta C+J)$,
$\text{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$, $\|CX_{C}-X\|^{2}_{F}$,
$\text{tr}(C^{T}BC)$, respectively. We prove this in Appendix B. After
ignoring the constant term, the majorised problem of equation 4 is

$\displaystyle\underset{C\in\mathcal{S}_{C}}{\min}$ $\displaystyle\
\frac{1}{2}C^{T}C-C^{T}\Big{(}C^{t}-\frac{1}{L}\nabla f(C^{t})\Big{)}$ (8)

Equation 8 is the majorized problem of equation 4. The optimal solution to
equation 8, found by using Karush–Kuhn–Tucker (KKT) optimality conditions is
(Proof is deferred to the Appendix B):

$\displaystyle C^{t+1}$ $\displaystyle=\Big{(}C^{t}-\frac{1}{L}\nabla
f\big{(}C^{t}\big{)}\Big{)}^{+}$ (9) $\displaystyle\text{where, }\nabla
f\big{(}C^{t}\big{)}$ $\displaystyle=-2\gamma\Theta C^{t}({C^{t}}^{T}\Theta
C^{t}+J)^{-1}+\alpha(C^{t}X_{C}-X)X_{C}^{T}+2\Theta
C^{t}X_{C}X_{C}^{T}-\frac{\beta}{e}BC^{t}$ (10)

Update rule of $X_{C}$. Treating $C$ fixed and $X_{c}$ as variable. The
subproblem for updating $X_{c}$ is

$\displaystyle\begin{array}[]{ll}\underset{\tilde{X}}{\min}f(\tilde{X})=\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})+\frac{\alpha}{2}\|X-CX_{C}\|_{F}^{2}\end{array}$ (12)

The closed form solution of problem equation 12 can be obtained by putting the
gradient of $f(\tilde{X})$ to zero.

$\displaystyle X_{C}^{t+1}=\Big{(}\frac{2}{\alpha}C^{T}\Theta
C+C^{T}C\Big{)}^{-1}C^{T}X$ (13)

Algorithm 1 Q-MAGC Algorithm

1:$G(X,\Theta),\alpha,\beta,\gamma,\lambda$

2:$t\leftarrow 0$

3:while Stopping Criteria not met do

4: Update $C^{t+1}$ as in Eqn. 9

5: Update $X_{C}^{t+1}$ as in Eqn. 13

6: $t\leftarrow t+1$

7:end while

8:return $C^{t},X_{C}^{t}$

Convergence Analysis.

######  Theorem 1.

The sequence $\\{C^{t+1},X_{C}^{t+1}\\}$ generated by Algorithm 1 converges to
the set of Karush–Kuhn–Tucker (KKT) optimality points for Problem 4

###### Proof.

The detailed proof can be found in the Appendix Section C. ∎

Complexity Analysis. The worst-case time complexity of a loop (i.e. one epoch)
in Algorithm 1 is $\mathcal{O}(p^{2}k+pkn)$ because of the matrix
multiplication in the update rule of $C$ equation 9. Here, $k$ is the number
of clusters and $n$ is the feature dimension. Note that $k$ is much smaller
than both $p$ and $n$. This makes our method much faster than previous
optimization based methods and faster than GCN-based clustering methods which
have complexities around $\mathcal{O}(p^{2}n+pn^{2})$. We discuss this more in
Appendix Section M.

Consistency Analysis on Degree Corrected Stochastic Block Models (DC-SBM).

A graph $G(V,E)$ is generated by a DC-SBM with $p$ nodes which belong to $k$
classes. Following the setup in Zhao et al. (2012), each node $v_{i}$ is
associated with a label-degree pair $(y_{i},t_{i})$ drawn from a discrete
joint distribution $\Pi_{k\times m}$. Here, $t_{i}$ takes values $0\leq
x_{1}\leq\dots\leq x_{m}$ with $\mathbb{E}[t_{i}]=1$, and $y_{i}\in[k]$.

Additionally, we have a symmetric $k\times k$ matrix $P$ that specifies the
probabilities of inter-cluster edges. The edges between each node pair
$(v_{i},v_{j})$ are sampled as independent Bernoulli trials with probability
$\mathbb{E}[A_{ij}]=t_{i}t_{j}P_{y_{i}y_{j}}$. To ensure
$\mathbb{E}[A_{ij}]<1$, DC-SBMs must satisfy $x_{m}^{2}(\max_{i,j}P_{ij})\leq
1$. The matrix $P$ is allowed to scale with the number of nodes $p$ and is
reparameterized as $P_{p}=\rho_{p}P$, where $\rho_{p}=\Pr[A_{ij}=1]\rightarrow
0$ as $p$ increases, and $P$ is fixed. The average degree of the DC-SBM graph
is defined as $\lambda_{p}=p\rho_{p}$.

Let $\hat{y}_{i}$ denote the predicted cluster for node $v_{i}$. The
assignment matrix $C$ is the one-hot encoding of $y$, such that
$C_{i}=\text{one-hot}(y_{i})$.

Next, we define $O(e)\in\mathbb{Z}_{\geq 0}^{k\times k}$ as the inter-cluster
edge count matrix for a given cluster assignment $e\in[k]^{p}$. Here,
$O_{ql}(e)=\sum_{ij}A_{ij}\mathbbm{1}\\{e_{i}=q,e_{j}=l\\}$ denotes the number
of edges between clusters $q$ and $l$. Additionally, we represent the class
frequency distribution as $\pi\in[0,1]^{k}$, where $\pi_{i}$ indicates the
fraction of nodes assigned to cluster $i$.

######  Definition 1.

(Strong and Weak Consistency). The clustering objective is defined to be
strongly consistent if $\lim\limits_{p\to\infty}Pr[\hat{y}=y]\to 1$. A weaker
notion of consistency is defined by
$\lim\limits_{p\to\infty}Pr\bigg{[}\frac{1}{p}\sum\limits_{i=1}^{p}\mathbbm{1}\\{\hat{y}\neq
y\\}<\epsilon\bigg{]}\to 1\ \ \forall\ \epsilon>0.$

Zhao et al. (2012) prove the consistency of multiple clustering objectives
under the DC-SBM, including Newman–Girvan modularity (Theorem 3.1). They also
provide a general theorem (Theorem 4.1) for checking consistency under DC-SBMs
for any criterion $\mathcal{L}(e)$, which can be expressed as
$\mathcal{L}(e)=F\left(\frac{O(e)}{\mu_{p}},\pi\right)$, where
$\mu_{p}=p^{2}\rho_{p}$ and $e\in[k]^{p}$ is a cluster assignment.

The consistency of an objective is evaluated by determining if the population
version of $\mathcal{L}(e)$ is maximized by the true cluster assignment $y$.
The population version of $\mathcal{L}(e)$ is obtained by taking its
conditional expectations given $y$ and $t$. We consider an array
$S\in\mathbb{R}^{k\times k\times m}$ and define a matrix
$H(S)\in\mathbb{R}^{k\times k}$ as
$H_{kl}(S)=\sum_{abuv}x_{u}x_{v}P_{ab}S_{kau}S_{lbv}$. Additionally, we define
a vector $h(S)\in\mathbb{R}^{k}$ as $h_{k}(S)=\sum_{au}S_{kau}$. Here, $H(S)$
and $h(S)$ denote the population versions of $O(e)$ and $\pi$, respectively,
such that $\frac{1}{\mu_{p}}\mathbb{E}[O|C,t]=H(S)$ and
$\mathbb{E}[\pi|C,t]=h(S)$.

######  Lemma 1.

(Theorem 4.1 from Zhao et al. (2012)) For any
$\mathcal{L}(e)=F(\frac{O(e)}{\mu_{p}},f(e))$, if $\pi,P,F$ is uniquely
maximized by $S=\mathbb{D}$ and satisfies the regularity conditions, then
$\mathcal{L}$ is strongly consistent under DC-SBMs
$\frac{\lambda_{p}}{logp}\rightarrow\infty$ and weakly consistent if
$\lambda_{p}\rightarrow\infty$:

  * 1.

$F$ is Lipschitz in its arguments ($H(S),h(S)$)

  * 2.

The directional derivative
$\frac{\partial^{2}}{\partial\varepsilon^{2}}F(M_{0}+\varepsilon(M_{1}-M_{0}),\mathbf{t_{0}}+\varepsilon(\mathbf{t_{1}}-\mathbf{t_{0}}))\big{|}_{\varepsilon=0^{+}}$
is continuous in $(M_{1},\mathbf{t_{1}})$ for all $(M_{0},\mathbf{t_{0}})$ in
a neighborhood of $(H(\mathbb{D}),\pi)$

  * 3.

With $G(S)=F(H(S),h(S))$, $\quad\frac{\partial
G((1-\varepsilon)\mathbb{D}+\varepsilon
S)}{\partial\varepsilon}|_{\varepsilon=0^{+}}<-C<0\ \forall\ \pi,P$

It is important to note that the paper focuses on maximizing an objective
function. Therefore, we consider the negative of our loss function,
$-\mathcal{L}_{MAGC}$. Additionally, we only consider the objective function
defined in Equation equation 4, as regularizers can be adjusted for different
downstream tasks.

######  Theorem 2.

Under the DC-SBM, $\mathcal{L}_{MAGC}$ is strongly consistent when
$\lambda_{p}/\log p\rightarrow\infty$ and weakly consistent when
$\lambda_{p}\rightarrow\infty$.

###### Proof.

The proof is divided into two parts. First, we demonstrate that
$\mathcal{L}_{MAGC}$ can be expressed in terms of $O$ and $\pi$. Second, we
verify that, in this form, $\mathcal{L}_{MAGC}$ satisfies the regularity
conditions specified in Theorem 4.1 of Zhao et al. (2012). It is important to
note that DC-SBMs do not consider node features, so we treat $X$ and $X_{C}$
as constants. We start by noting that $O=C^{T}AC$, i.e., the edge count matrix
$O$ is the class-transformed version of the adjacency matrix $A$. This has a
good intuition behind it, as $C$ is the mapping from nodes to classes, and it
can also be proved easily.

$\\{C^{T}AC\\}_{ql}=(C^{T})_{q}AC_{:l}=\sum_{i=1}^{p}\sum_{j=1}^{p}(C^{T})_{qi}A_{ij}C_{jl}=\sum_{i=1}^{p}\sum_{j=1}^{p}C_{iq}A_{ij}C_{jl}$
(14)

$C_{iq}$ equals 1 only when node $i$ is in cluster $q$, and $C_{jl}$ equals 1
only when node $j$ is in cluster $l$. Therefore, this expression simplifies to
the definition of $O$.

$\\{C^{T}AC\\}_{ql}=\sum_{i=1}^{p}\sum_{j=1}^{p}A_{ij}\mathbbm{1}\\{y_{i}=q,y_{j}=l\\}$
(15)

This can be viewed as the _cluster-level adjacency matrix_ , which quantifies
the edges between clusters. Next, we formalize the relationship between $O$
and $C^{T}DC$, where $D=\text{diag}(d)$ is the degree matrix.

$\\{C^{T}DC\\}_{ql}=\sum_{i=1}^{p}\sum_{j=1}^{p}C_{iq}D_{ij}C_{jl}=\sum_{i=1}^{p}C_{iq}d_{i}C_{il}\
\ \ [D_{ij}=0\text{ if }i\neq j]\\\ $ (16)

Since the rows of $C$ are orthonormal, $C_{iq}C_{il}=1$ if and only if $q=l$.
Thus, we obtain the following expression for the diagonal entries of
$C^{T}DC$:

$\\{C^{T}DC\\}_{qq}=\sum_{i=1}^{p}C_{iq}^{2}d_{i}=\sum_{i=1}^{p}C_{iq}d_{i}\ \
\ \bigg{[}\because C_{ij}\in\\{0,1\\}\bigg{]}\\\ $ (17)

Now, we consider $\sum_{l=1}^{k}\\{C^{T}AC\\}_{ql}$

$=\sum_{i=1}^{p}\sum_{j=1}^{p}\bigg{(}C_{iq}A_{ij}\sum_{l=1}^{k}C_{jl}\bigg{)}=\sum_{i=1}^{p}\sum_{j=1}^{p}C_{iq}A_{ij}\cdot
1=\sum_{i=1}^{p}C_{iq}d_{i}$ (18)

Thus, we showed that $C^{T}DC=\text{diag}(\sum_{l=1}^{k}[C^{T}AC]_{k})$, which
intuitively represents the _cluster-level degree matrix_ because it is the
(diagonal of) row/column-sum of the $C^{T}AC$ or cluster-level adjacency
matrix. We now express the Laplacian of the coarsened graph as a function of
$O$, $\Theta_{C}(O)=C^{T}\Theta
C=C^{T}(D-A)C=\text{diag}(\sum_{l=1}^{k}O_{l})-O$. This implies that the terms
$\textrm{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$, $\textrm{tr}(C^{T}BC)$, and
$\log\det(C^{T}\Theta C+J)$ in our clustering objective are functions of
$\Theta_{C}(O)$. Using Theorem 4.1 from Zhao et al. (2012), we conclude that
$\mathcal{L}_{\text{MAGC}}$ must be uniquely minimized at any point
$(y^{*},A^{*})$ s.t. $\mathbb{E}_{p}[\pi(y_{p})]=\pi(y^{*})$ and
$\mathbb{E}_{p}[A^{(p)}]=A^{*}$. The Laplacian $\Theta$ can be written as a
function of the adjacency $\Theta_{ij}=\sum_{j}A^{*}_{ij}-A^{*}_{ij}$
(assuming no self-loops).

For the modularity term, this has been proven in Theorem 3.1 of Zhao et al.
(2012). For the constraint relaxation term $\frac{1}{2}\|CX_{C}-X\|_{F}^{2}$,
we see that $X_{C}$ is a function of only $C(y^{*})$ and $\Theta(A^{*})$ from
it’s update rule equation 13, so we get
$\text{term}_{\alpha}(y^{*},\Theta)=\frac{1}{2}\|C(y^{*})X_{C}(y^{*},A^{*})-X\|_{F}^{2}$.
$\text{tr}(X_{C}C^{T}\Theta CX_{C})$ is already just a function of
$X_{C},C,\Theta$ which we have shown above to be functions of only
$y^{*},A^{*}$. By the same logic, $-\log\det(C^{T}\Theta C+J)$ is also a
function of only $y^{*},A^{*}$ ($J$ is a constant =
$\frac{1}{k}\mathbf{1}_{k\times k}$). Thus, $\mathcal{L}_{\text{MAGC}}$ can be
written as required above.

Next, we compute the population version of $\mathcal{L}_{MAGC}$.

$\displaystyle F(H(S),h(S))$
$\displaystyle=\mathbb{E}[-\mathcal{L}_{MAGC}|C,t]=\underbrace{f_{1}(H(S),h(S))}_{\text{smoothness}}+\underbrace{f_{2}(H(S),h(S))}_{\text{log-
determinant}}+\underbrace{f_{3}(H(S),h(S))}_{\text{modularity}}$ (19)
$\displaystyle\text{where, }f_{1}(H(S),h(S))$
$\displaystyle=\mu_{p}\text{tr}(X_{C}^{T}[H(S)-\text{diag}(\sum_{j=1}^{k}H(S)_{ij})]X_{C})$
(20) $\displaystyle f_{2}(H(S),h(S))$
$\displaystyle=\frac{\mu_{p}}{b}\text{tr}\bigg{(}H(S)-\text{diag}(\sum_{j=1}^{k}H(S)_{ij})\bigg{)}-\frac{1}{b}+k-k\log
b$ (21) $\displaystyle f_{3}(H(S),h(S))$
$\displaystyle=\frac{\text{tr}(H(S))}{\tilde{P_{0}}}-\frac{\sum_{i=1}^{k}(\sum_{j=1}^{k}H(S)_{ij})^{2}}{\tilde{P_{0}}^{2}}$
(22) $\displaystyle\text{and, }\tilde{P_{0}}$
$\displaystyle=\sum_{ab}\tilde{\pi}_{a}\tilde{\pi}_{b}P_{ab}\quad;\quad\tilde{\pi}_{a}=\sum_{u}x_{u}\Pi_{au}\text{
with }\sum_{a}\tilde{\pi}_{a}=1,\text{ since }\mathbb{E}[t_{i}]=1.$ (23)

For brevity, we denote $F(H(S),h(S))$ by $F(S)$. Now for the second part of
the proof, it is easy to prove that $F(S)$ satisfies the regularity conditions
stated in Lemma 1. The detailed proof is deferred to Appendix D. This
concludes the proof that $\mathcal{L}_{MAGC}$ is strongly and weakly
consistent under the DC-SBM. ∎

We also demonstrate experimental consistency of the optimization objective
equation 4 in Appendix I, resulting in complete recovery of the labels. This
consistency extends to the objective with additional regularization terms, as
detailed therein.

##  5 Integrating with GNNs

Algorithm 2 Q-GCN Algorithm

1:$G(X,\Theta),\alpha,\beta,\gamma,\lambda$

2:$k$ GCN layers - weights $W^{(1..k)}$

3:$t\leftarrow 0$

4:while Stopping Criteria not met do

5: $C^{t+1}\leftarrow\text{GCN}(X,A)$

6: $X_{C}^{t+1}\leftarrow C^{\dagger}X$

7: $t\leftarrow t+1$

8: $\mathcal{L}\leftarrow\mathcal{L}_{\text{MAGC}}(C,X)$

9:end while

10:return $C^{t},X_{C}^{t}$

Our optimization framework integrates seamlessly with Graph Neural Networks
(GNNs) by incorporating the optimization objective equation 4 into the loss
function. This integration can be minimized using gradient descent. We
demonstrate the effectiveness of this approach on several popular GNN
architectures, including Graph Convolutional Networks (GCNs) (Kipf & Welling,
2016b), Variational Graph Auto-Encoders (VGAEs) (Kipf & Welling, 2016a), and a
variant known as Gaussian Mixture Model VGAE (GMM-VGAE) (Hui et al., 2020).

We iteratively learn the matrix $C$ using gradient descent. From equation 1,
we update $X_{C}$ using the relation $X_{C}=C^{\dagger}X$. It is important to
note that the loss term $CX_{C}-X$ will not necessarily be zero. This arises
because we use a "soft" version of $C$ in the loss function to enable gradient
flow, while a "hard" version of $C$ is used in the update step. This method
ensures that $C$ naturally becomes harder and exhibits a higher prediction
probability.

Q-GCN. We integrate our loss function ($\mathcal{L}_{MAGC}$, as defined in
equation 4) into a simple three-layer Graph Convolutional Network (GCN) model.
The soft cluster assignments $C$ are learned as the output of the final GCN
layer. The detailed architecture and loss function are illustrated in Figure
1.

Q-VGAE. The VGAE loss can be written as
$\mathcal{L}_{VGAE}=\lambda_{recon}\underbrace{\mathbb{E}_{q(Z|X,A)}[\log
p(\hat{A}|Z)]}_{\text{Reconstruction
Error}}-\lambda_{kl}\underbrace{\textrm{KL}[q(Z|X,A)\ ||\
p(Z)]}_{\text{Kullback-Leibler divergence}}$

Algorithm 3 Q-VGAE/Q-GMM-VGAE Algorithm
$Z$ is the latent space of the VGAE and $\hat{A}$ is the reconstructed
adjacency matrix.

1:$G(X,\Theta),\alpha,\beta,\gamma,\lambda$

2:Variational Encoder - VarEnc (GCN, $\mu$, $\sigma$)

3:Prediction Head - Pred (GCN or GMM)

4:$t\leftarrow 0$

5:while Stopping Criteria not met do

6: $Z\leftarrow\text{VarEnc}_{\mu,\sigma}(X,A)$

7: $\hat{A}\leftarrow ZZ^{T}$

8: $C^{t+1}\leftarrow\text{Pred}(Z,A)$

9: $X_{C}^{t+1}\leftarrow C^{\dagger}X$

10: $t\leftarrow t+1$

11:
$\mathcal{L}\leftarrow\mathcal{L}_{\text{MAGC}}(C,X)+\mathcal{L}_{\text{VGAE}}(X,A,Z,\hat{A})$

12:end while

13:return $C^{t},X_{C}^{t}$

where, $Z$ represents the latent space of the VGAE, $\hat{A}$ is the
reconstructed adjacency matrix, and $\lambda_{recon}$ and $\lambda_{kl}$ are
hyperparameters. We add a GCN layer on top of this architecture, which takes
$Z$ as input and predicts $C$. A detailed summary of the VGAE loss terms is
provided in Appendix Section G. For the VGAE, we minimize the sum of three
losses: the reconstruction loss, the KL-divergence loss, and our loss. This
combined loss function is expressed as:
$\mathcal{L}_{Q-VGAE}=\mathcal{L}_{MAGC}+\mathcal{L}_{VGAE}$

Q-GMM-VGAE. This variant of VGAE incorporates a Gaussian Mixture Model (GMM)
in the latent space to better capture data distributions. This approach is
effective because it minimizes the evidence lower bound (ELBO) or variational
lower bound (Hui et al., 2020; Kingma & Welling, 2014; Kipf & Welling, 2016a)
using multiple priors, rather than a single Gaussian prior as in standard
VGAE. Hui et al. (2020) use a number of priors equal to the number of
clusters.

Figure 1: a) Architecture of Q-GCN. We want to train the encoder to learn the
soft cluster assignment matrix $C$. The coarsened features $X_{C}$ are
obtained using the relation $X_{C}^{t+1}={C^{t+1}}^{\dagger}X$. Finally, our
proposed MAGC loss is then computed using $C$ and $X_{C}$.
b) Architecture of Q-VGAE/Q-GMM-VGAE. The three-layer GCN encoder takes $X$
and $A$ as inputs to learn the latent representation $Z$ of the graph. $Z$ is
then passed through an inner-product decoder to reconstruct the adjacency
matrix $\hat{A}$. The reconstruction loss is calculated between $\hat{A}$ and
$A$, and the KL-divergence is applied to $Z$. In Q-VGAE (or Q-GMM-VGAE), $Z$
is also passed through a GCN layer (or GMM) to output the soft cluster
assignments $C$. The MAGC loss is then computed in the same manner as in
Q-GCN.

##  6 Experiments

###  6.1 Benchmark Datasets and Baselines

We evaluate our method on a diverse set of datasets, including small
attributed datasets like Cora and CiteSeer, larger datasets like PubMed, and
unattributed datasets such as Airports (Brazil, Europe, and USA).
Additionally, we test our method on very large datasets like
CoauthorCS/Physics, AmazonPhoto/PC, and ogbn-arxiv. A detailed summary of all
the datasets used is provided in Appendix E.

To assess the performance of our method, we compare it against three types of
state-of-the-art methods based on the input and architecture type: methods
that use only node attributes, methods that use only graph structure, and
methods that use both graph structure and node attributes. The last category
is further subdivided into graph coarsening methods, GCN-based architectures,
VGAE-based architectures and contrastive methods, and heavily modified VGAE
architectures. This comprehensive evaluation allows us to demonstrate the
robustness and versatility of our approach across various data and model
configurations.

###  6.2 Metrics

To evaluate the performance of our method, we utilize label alignment metrics
that compare ground truth node labels with cluster assignments. Specifically,
we measure Normalised Mutual Information (NMI), Adjusted Rand Index (ARI), and
Accuracy (ACC), with higher values indicating superior performance. For
detailed explanations of these metrics, please refer to Appendix E.

We selected Normalized Mutual Information (NMI) as the primary metric for
evaluating model performance based on its prominence in graph clustering
literature and its comprehensive ability to assess the quality of our cluster
assignments.

Training Details. Training details are available in the Appendix F.

###  6.3 Attributed Graph Clustering

We present our key results on the real datasets Cora, CiteSeer, and PubMed in
Table 1. Our proposed method outperforms all existing methods in terms of NMI
and demonstrates competitive performance in Accuracy and ARI. The best models
were selected based on NMI scores. Results for very large datasets are
provided in Appendix J. Unlike some methods such as S3GC (Devvrit et al.,
2022), which use randomly-sampled batches that can introduce bias by breaking
community structure, we perform full-batch training by passing the entire
graph. For extremely large graphs, such as ogbn-arxiv, we have also utilized
batching.

Cora CiteSeer PubMed Method ACC $\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC
$\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC $\uparrow$ NMI $\uparrow$ ARI
$\uparrow$ K-means 34.7 16.7 25.4 38.5 17.1 30.5 57.3 29.1 57.4 Spectral
Clustering 34.2 19.5 30.2 25.9 11.8 29.5 39.7 3.5 52.0 DeepWalk (Perozzi et
al., 2014a) 46.7 31.8 38.1 36.2 9.7 26.7 61.9 16.7 47.1 Louvain (Blondel et
al., 2008b) 52.4 42.7 24.0 49.9 24.7 9.2 30.4 20.0 10.3 GAE [NeurIPS’16] (Kipf
& Welling, 2016a) 61.3 44.4 38.1 48.2 22.7 19.2 63.2 24.9 24.6 DGI [ICLR’19]
(Veličković et al., 2019) 71.3 56.4 51.1 68.8 44.4 45.0 53.3 18.1 16.6 GIC
[PAKDD’21] (Mavromatis & Karypis, 2021) 72.5 53.7 50.8 69.6 45.3 46.5 67.3
31.9 29.1 DAEGC [IJCAI’19] (Wang et al., 2019) 70.4 52.8 49.6 67.2 39.7 41.0
67.1 26.6 27.8 GALA [ICCV’19] (Park et al., 2019) 74.5 57.6 53.1 69.3 44.1
44.6 69.3 32.1 32.1 AGE [KDD’20] (Cui et al., 2020) 73.5 57.5 50.0 69.7 44.9
34.1 71.1 31.6 33.4 DCRN [AAAI’22] (Liu et al., 2022) 61.9 45.1 33.1 70.8 45.8
47.6 69.8 32.2 31.4 FGC [JMLR’23] (Kumar et al., 2023) 53.8 23.2 20.5 54.2
31.1 28.2 67.1 26.6 27.8 Q-MAGC (Ours) 65.8 51.8 42.0 65.9 40.8 40.1 66.7 32.8
27.9 Q-GCN (Ours) 71.6 58.3 53.6 71.5 47.0 49.1 64.1 32.1 26.5 SCGC [IEEE
TNNLS’23] (Liu et al., 2023a) 73.8 56.1 51.7 71.0 45.2 46.2 - - - MVGRL
[ICML’20] (Hassani & Khasahmadi, 2020) 73.2 56.2 51.9 68.1 43.2 43.4 69.3 34.4
32.3 VGAE [NeurIPS’16] (Kipf & Welling, 2016a) 64.7 43.4 37.5 51.9 24.9 23.8
69.6 28.6 31.7 ARGA [IJCAI’18] (Pan et al., 2018) 64.0 35.2 61.9 57.3 34.1
54.6 59.1 23.2 29.1 ARVGA [IJCAI’18] (Pan et al., 2018) 63.8 37.4 62.7 54.4
24.5 52.9 58.2 20.6 22.5 R-VGAE [IEEE TKDE’22] (Mrabah et al., 2022) 71.3 49.8
48.0 44.9 19.9 12.5 69.2 30.3 30.9 Q-VGAE (Ours) 72.7 58.6 49.6 66.1 47.4 50.2
64.3 32.6 28.0 VGAECD-OPT [Entropy’20] (Choong et al., 2020) 27.2 37.3 22.0
51.8 25.1 15.5 32.2 25.0 26.1 Mod-Aware VGAE [NN’22] (Salha-Galvan et al.,
2022) 67.1 52.4 44.8 51.8 25.1 15.5 - 30.0 29.1 GMM-VGAE [AAAI’20] (Hui et
al., 2020) 71.9 53.3 48.2 67.5 40.7 42.4 71.1 29.9 33.0 R-GMM-VGAE [IEEE
TKDE’22] (Mrabah et al., 2022) 76.7 57.3 57.9 68.9 42.0 43.9 74.0 33.4 37.9
Q-GMM-VGAE (Ours) 76.2 58.7 56.3 72.7 47.4 48.8 69.0 34.8 34.0

Table 1: Comparison of all methods on attributed datasets. We classify the
baselines into three primary groups: the first includes traditional clustering
algorithms and GNN-free methods, serving as relevant baselines for our Q-MAGC
method. The second category compares GNN-based methods with Q-GCN, and the
last, comprising the most performant methods, consists of VGAE-based baselines
for Q-VGAE and Q-GMM-VGAE.

###  6.4 Non-Attributed Graph Clustering

For non-attributed graphs, we use a one-hot encoding of the degree vector as
features. While this is a basic approach compared to learning-based methods
like DeepWalk (Perozzi et al., 2014b) and node2vec (Grover & Leskovec, 2016),
it ensures a fair comparison since other methods also use this feature
representation. Our results in Table 2(a) demonstrate that our method achieves
competitive or superior performance in terms of NMI, even for non-attributed
datasets.

###  6.5 Ablation Studies

Comparison of running times In Fig 2(a), we compare the running times of our
method with other baselines. Our method consistently takes less than half the
time across all datasets. Notably, on PubMed (large dataset), state-of-the-art
methods GMM-VGAE and R-GMM-VGAE (unmodified) require approximately 60 minutes
for clustering, whereas our Q-GMM-VGAE delivers superior performance in under
15 minutes, representing a $75\%$ reduction. Additionally, Q-MAGC runs even
faster, completing in just 6 minutes on PubMed, and achieves approximately
$90\%$ of the performance.

Brazil Europe USA Method ACC $\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC
$\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC $\uparrow$ NMI $\uparrow$ ARI
$\uparrow$ GAE [NeurIPS’16] 62.6 37.8 30.8 47.6 19.9 12.7 43.9 13.6 11.8 DGI
[ICLR’19] 64.9 31.0 30.4 48.6 16.1 12.3 52.2 22.9 21.7 GIC [PAKDD’21] 40.5
23.5 14.1 40.4 9.4 6.2 49.7 22.1 19.9 DAEGC [AAAI’19] 71.0 47.4 41.2 53.6 30.9
23.3 46.4 27.2 18.4 Q-GCN (Ours) 51.1 31.9 23.7 45.5 30.8 25.1 43.8 19.1 14.8
VGAE [NeurIPS’16] 64.1 38.0 30.7 49.9 23.5 16.7 45.8 23.1 15.7 Q-VGAE (Ours)
50.1 35.0 19.8 46.6 19.5 17.5 46.2 19.5 16.9 GMM-VGAE [AAAI’20] 70.2 46.0 41.9
53.1 31.1 24.4 48.1 21.9 13.2 R-GMM-VGAE [IEEE TKDE’22] 73.3 45.6 42.5 57.4
31.4 25.8 50.8 23.1 15.3 Q-GMM-VGAE (Ours) 68.4 46.0 42.4 47.9 32.2 23.5 46.6
23.1 13.1

(a) Comparison of all methods on non-attributed datasets using degree. (a)
Comparison of running times of methods.
Note that the scale for PubMed is in minutes (right axis), whereas for Cora
and CiteSeer is in seconds.

Modularity Metric Comparison We treat modularity as a metric and measure the
gains observed in modularity over other baselines on the Cora, CiteSeer, and
PubMed datasets. We report two types of graph-based metrics: modularity
$\mathcal{Q}$ and conductance $\mathcal{C}$, which do not require labels.
Conductance measures the fraction of total edge volume pointing outside the
cluster, with $\mathcal{C}$ being the average conductance across all clusters,
where a lower value is preferred.

From Table 3(a), we observe that although DMoN (Tsitsulin et al., 2023)
achieves the highest modularity, our method attains significantly higher NMI.
For CiteSeer, we achieve a $40\%$ improvement in NMI with only an $8\%$
decrease in modularity, positioning us closer to the ground truth.
Additionally, our methods outperform their foundational counterparts, with
Q-MAGC outperforming FGC, and Q-VGAE outperforming VGAE. Although modularity
is a valuable metric to optimize, the maximum modularity labeling of a graph
does not always correspond to the ground truth labeling. Therefore, it is
crucial to include the other terms in our formulation. While optimizing
modularity helps approach the optimal model parameters (where NMI would be 1),
it can deviate slightly. The additional terms in our formulation correct this
trajectory, ensuring more accurate results.

Figure 3: Evolution of the latent space of a) Q-VGAE and b) Q-GMM-VGAE over
time for Cora. Colors represent cluster assignments.

Cora CiteSeer PubMed $\mathcal{C}\downarrow$ $\mathcal{Q}\uparrow$ NMI
$\uparrow$ $\mathcal{C}\downarrow$ $\mathcal{Q}\uparrow$ NMI $\uparrow$
$\mathcal{C}\downarrow$ $\mathcal{Q}\uparrow$ NMI $\uparrow$ DMoN 12.2 76.5
48.8 5.1 79.3 33.7 17.7 65.4 29.8 FGC 58.4 25 23.1 41.6 41.1 31 21.6 44.1 20.5
Q-MAGC 13.3 72.5 51.7 16.8 64.9 40.16 26 40.3 28.1 Q-GCN 13.6 73.3 58.3 5.8
74.5 46.7 8.27 55 31.5 VGAE 17.6 60.8 38.1 12.8 55.8 21 13.5 45.8 26.9 Q-VGAE
9.5 71.5 58.4 4.6 72.4 47.3 9.4 52.12 31.8

(a) Comparison of modularity and conductance at the best NMI with DMoN. Note
that DMoN is optimizing only modularity, whereas we are optimizing other
important terms as well, as mentioned in Eqn 4, and thus gain a lot on NMI by
giving up a small amount of modularity, making us closer to the ground truth.
(a) Impact of active parameters on clustering performance. All terms in the
loss equation 4 are represented by their parameters, for example the
modularity term is represented by $\beta$. Additionally, $\omega$ represents
the non-parameterized term.

Importance of and Evolution of different loss terms We analyze the evolution
of the different loss terms during training, and also try to measure the
impact of each term separately by removing terms from the loss one by one, as
shown in Appendix N. Also, we found that $||CX_{C}-X||_{F}^{2}$ is the most
sensitive to change in its weight $\alpha$, followed by the terms related to
$\gamma$, $\beta$ and then $\lambda$. This makes sense because if that
constraint(relaxation) is not being met, then $C$ would have errors. Even
though some of the terms do the heavy lifting, the other regularization terms
do contribute to performance and more importantly, change the nature of $C$ :
The smoothness term corresponds to smoothness of signals in the graph being
transferred to the coarsened graph which encourages local "patches"/groups of
$C$ to belong to the same cluster. The term $\gamma$ ensures that the
coarsened graph is connected - i.e. preserving inter-cluster relations, which
simple contrastive methods destroy; this affects $C$ by making it so that
$\Theta_{C}$ has minimal multiplicity of 0-eigenvalues.

Visualization of the latent space In Figure 3, we visualize how the latent
space of the Q-VGAE and Q-GMM-VGAE changes over time for the Cora dataset.
Plots for the rest of the datasets can be found in the Appendix H. We use UMAP
(Uniform Manifold Approximation and Projection) (McInnes et al., 2018) for
dimensionality reduction.

##  7 Discussion

Performance of. Our methods consistently outperform their counterparts. Q-MAGC
demonstrates significant superiority over traditional algorithms and GNN-free
approaches. These can be attributed to the inclusion of modularity in our
objective, as well as the smoothness and connectedness terms, making Q-MAGC
more effective than previous methods.

To compete with GNN-based methods, we leverage their powerful message-passing
mechanisms. Q-GCN surpasses other GCN/simple GNN-based methods, such as AGE
and DCRN, not only because it optimizes modularity but also due to its focus
on smoothness and inter-cluster relationships in the coarsened graph.
Similarly, Q-VGAE achieves better results than other VGAE/complex GNN-based
methods for the same reasons.

Efficiency. Our methods are efficient in both complexity and implementation,
as demonstrated in the Complexity Analysis 1 and the Ablation - Comparison of
running times.

Limitations. If the modularity of a graph, calculated from the ground truth
labeling, is low, our method may not perform optimally and can be slower to
converge since maximizing modularity in such cases is not ideal. However, our
method still manages to match or surpass state-of-the-art methods on the
Airports dataset, where all graphs exhibit low modularity based on ground
truth labels. This is due to the other terms in our weighted loss function,
which become the primary contributors to performance.

## References

  * Arvanitidis et al. (2018) Georgios Arvanitidis, Lars Kai Hansen, and Søren Hauberg.  Latent space oddity: on the curvature of deep generative models.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=SJzRZ-WCZ.
  * Bhowmick et al. (2023) Aritra Bhowmick, Mert Kosan, Zexi Huang, Ambuj Singh, and Sourav Medya.  Dgcluster: A neural framework for attributed graph clustering via modularity maximization, 2023.
  * Bianchi et al. (2020) Filippo Maria Bianchi, Daniele Grattarola, and Cesare Alippi.  Spectral clustering with graph neural networks for graph pooling.  In _Proceedings of the 37th International Conference on Machine Learning_ , ICML’20. JMLR.org, 2020.
  * Bickel & Chen (2009) Peter J. Bickel and Aiyou Chen.  A nonparametric view of network models and newman–girvan and other modularities.  _Proceedings of the National Academy of Sciences_ , 106(50):21068–21073, 2009.  doi: 10.1073/pnas.0907096106.  URL https://www.pnas.org/doi/abs/10.1073/pnas.0907096106.
  * Blondel et al. (2008a) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre.  Fast unfolding of communities in large networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2008(10):P10008, oct 2008a.  doi: 10.1088/1742-5468/2008/10/P10008.  URL https://dx.doi.org/10.1088/1742-5468/2008/10/P10008.
  * Blondel et al. (2008b) Vincent D Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre.  Fast unfolding of communities in large networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2008(10):P10008, October 2008b.  ISSN 1742-5468.  doi: 10.1088/1742-5468/2008/10/p10008.  URL http://dx.doi.org/10.1088/1742-5468/2008/10/P10008.
  * Bogdanov & Shchur (2021) Igor Bogdanov and Vladimir Shchur.  Variational autoencoders with euclidean and hyperbolic latent spaces for population genetics.  In _2021 XVII International Symposium "Problems of Redundancy in Information and Control Systems" (REDUNDANCY)_ , pp. 91–94, 2021.  doi: 10.1109/REDUNDANCY52534.2021.9606448.
  * Brandes et al. (2008) Ulrik Brandes, Daniel Delling, Marco Gaertler, R. Gorke, Martin Hoefer, Zoran Nikoloski, and Dorothea Wagner.  On modularity clustering.  _Knowledge and Data Engineering, IEEE Transactions on_ , 20:172 – 188, 03 2008.  doi: 10.1109/TKDE.2007.190689.
  * Buterez et al. (2021) David Buterez, Ioana Bica, Ifrah Tariq, Helena Andrés-Terré, and Pietro Liò.  CellVGAE: an unsupervised scRNA-seq analysis workflow with graph attention networks.  _Bioinformatics_ , 38(5):1277–1286, 12 2021.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btab804.  URL https://doi.org/10.1093/bioinformatics/btab804.
  * Caron et al. (2018) Mathilde Caron, Piotr Bojanowski, Armand Joulin, and Matthijs Douze.  Deep clustering for unsupervised learning of visual features.  In _European Conference on Computer Vision_ , 2018.
  * Chen et al. (2020) Nutan Chen, Alexej Klushyn, Francesco Ferroni, Justin Bayer, and Patrick Van Der Smagt.  Learning flat latent manifolds with VAEs.  In Hal Daumé III and Aarti Singh (eds.), _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pp. 1587–1596. PMLR, 13–18 Jul 2020.  URL https://proceedings.mlr.press/v119/chen20i.html.
  * Cheng & Ma (2022) Yi Cheng and Xiuli Ma.  scGAC: a graph attentional architecture for clustering single-cell RNA-seq data.  _Bioinformatics_ , 38(8):2187–2193, 02 2022.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btac099.  URL https://doi.org/10.1093/bioinformatics/btac099.
  * Choong et al. (2020) Jun Jin Choong, Xin Liu, and Tsuyoshi Murata.  Optimizing variational graph autoencoder for community detection with dual optimization.  _Entropy_ , 22(2), 2020.  ISSN 1099-4300.  URL https://www.mdpi.com/1099-4300/22/2/197.
  * Connor et al. (2021) Marissa Connor, Gregory Canal, and Christopher Rozell.  Variational autoencoder with learned latent structure.  In Arindam Banerjee and Kenji Fukumizu (eds.), _Proceedings of The 24th International Conference on Artificial Intelligence and Statistics_ , volume 130 of _Proceedings of Machine Learning Research_ , pp. 2359–2367. PMLR, 13–15 Apr 2021.  URL https://proceedings.mlr.press/v130/connor21a.html.
  * Cui et al. (2020) Ganqu Cui, Jie Zhou, Cheng Yang, and Zhiyuan Liu.  Adaptive graph encoder for attributed graph embedding.  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, KDD ’20, pp. 976–985, New York, NY, USA, 2020. Association for Computing Machinery.  ISBN 9781450379984.  doi: 10.1145/3394486.3403140.  URL https://doi.org/10.1145/3394486.3403140.
  * Devvrit et al. (2022) Fnu Devvrit, Aditya Sinha, Inderjit S Dhillon, and Prateek Jain.  S3GC: Scalable self-supervised graph clustering.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho (eds.), _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=ldl2V3vLZ5.
  * Dhillon et al. (2007) Inderjit S. Dhillon, Yuqiang Guan, and Brian Kulis.  Weighted graph cuts without eigenvectors a multilevel approach.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 29(11):1944–1957, 2007.  doi: 10.1109/TPAMI.2007.1115.
  * Fey & Lenssen (2019) Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with pytorch geometric, 2019.  URL http://arxiv.org/abs/1903.02428.  cite arxiv:1903.02428.
  * Fiedler (1973) Miroslav Fiedler.  Algebraic connectivity of graphs.  _Czechoslovak Mathematical Journal_ , 23(2):298–305, 1973.  URL http://eudml.org/doc/12723.
  * Fortunato & Barthélemy (2007) Santo Fortunato and Marc Barthélemy.  Resolution limit in community detection.  _Proceedings of the National Academy of Sciences_ , 104(1):36–41, 2007.  doi: 10.1073/pnas.0605965104.  URL https://www.pnas.org/doi/abs/10.1073/pnas.0605965104.
  * Gao & Ji (2019) Hongyang Gao and Shuiwang Ji.  Graph u-nets.  In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 2083–2092. PMLR, 09–15 Jun 2019.  URL https://proceedings.mlr.press/v97/gao19a.html.
  * Grover & Leskovec (2016) Aditya Grover and Jure Leskovec.  Node2vec: Scalable feature learning for networks.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 855–864, New York, NY, USA, 2016. Association for Computing Machinery.  ISBN 9781450342322.  doi: 10.1145/2939672.2939754.  URL https://doi.org/10.1145/2939672.2939754.
  * Guimerà & Nunes Amaral (2005) Roger Guimerà and Luís A Nunes Amaral.  Functional cartography of complex metabolic networks.  _Nature_ , 433(7028):895–900, February 2005.  ISSN 0028-0836.  doi: 10.1038/nature03288.
  * Guimerà & Amaral (2005) Roger Guimerà and Luís Amaral.  Cartography of complex networks: Modules and universal roles.  _Journal of statistical mechanics (Online)_ , 2005:nihpa35573, 03 2005.  doi: 10.1088/1742-5468/2005/02/P02001.
  * Hamaguchi et al. (2017) Takuo Hamaguchi, Hidekazu Oiwa, Masashi Shimbo, and Yuji Matsumoto.  Knowledge transfer for out-of-knowledge-base entities : A graph neural network approach.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17_ , pp. 1802–1808, 2017.  doi: 10.24963/ijcai.2017/250.  URL https://doi.org/10.24963/ijcai.2017/250.
  * Hassani & Khasahmadi (2020) Kaveh Hassani and Amir Hosein Khasahmadi.  Contrastive multi-view representation learning on graphs.  In _Proceedings of International Conference on Machine Learning_ , pp. 3451–3461. 2020.
  * Hui et al. (2020) Binyuan Hui, Pengfei Zhu, and Qinghua Hu.  Collaborative graph convolutional networks: Unsupervised learning meets semi-supervised learning.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34(04):4215–4222, Apr. 2020.  doi: 10.1609/aaai.v34i04.5843.  URL https://ojs.aaai.org/index.php/AAAI/article/view/5843.
  * Karrer & Newman (2011) Brian Karrer and M. E. J. Newman.  Stochastic blockmodels and community structure in networks.  _Phys. Rev. E_ , 83:016107, Jan 2011.  doi: 10.1103/PhysRevE.83.016107.  URL https://link.aps.org/doi/10.1103/PhysRevE.83.016107.
  * Kingma & Welling (2014) Diederik P. Kingma and Max Welling.  Auto-Encoding Variational Bayes.  In _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , 2014.
  * Kipf & Welling (2016a) Thomas N. Kipf and Max Welling.  Variational graph auto-encoders, 2016a.
  * Kipf & Welling (2016b) Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  arXiv, 2016b.  doi: 10.48550/ARXIV.1609.02907.  URL https://arxiv.org/abs/1609.02907.
  * Kumar et al. (2023) Manoj Kumar, Anurag Sharma, and Sandeep Kumar.  A unified framework for optimization-based graph coarsening.  _Journal of Machine Learning Research_ , 24(118):1–50, 2023.  URL http://jmlr.org/papers/v24/22-1085.html.
  * Kumar et al. (2020) Sandeep Kumar, Jiaxi Ying, José Vinícius de M Cardoso, and Daniel P Palomar.  A unified framework for structured graph learning via spectral constraints.  _Journal of Machine Learning Research_ , 21(22):1–60, 2020.
  * Lee et al. (2019) Junhyun Lee, Inyeop Lee, and Jaewoo Kang.  Self-attention graph pooling.  In _Proceedings of the 36th International Conference on Machine Learning_ , 09–15 Jun 2019.
  * Liu et al. (2022) Yue Liu, Wenxuan Tu, Sihang Zhou, Xinwang Liu, Linxuan Song, Xihong Yang, and En Zhu.  Deep graph clustering via dual correlation reduction.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pp. 7603–7611, 2022.
  * Liu et al. (2023a) Yue Liu, Xihong Yang, Sihang Zhou, Xinwang Liu, Siwei Wang, Ke Liang, Wenxuan Tu, and Liang Li.  Simple contrastive graph clustering.  _IEEE Transactions on Neural Networks and Learning Systems_ , pp. 1–12, 2023a.  doi: 10.1109/TNNLS.2023.3271871.
  * Liu et al. (2023b) Yue Liu, Xihong Yang, Sihang Zhou, Xinwang Liu, Zhen Wang, Ke Liang, Wenxuan Tu, Liang Li, Jingcan Duan, and Cancan Chen.  Hard sample aware network for contrastive deep graph clustering.  In _Proc. of AAAI_ , 2023b.
  * Loukas (2019) Andreas Loukas.  Graph reduction with spectral and cut guarantees.  _Journal of Machine Learning Research_ , 20(116):1–42, 2019.  URL http://jmlr.org/papers/v20/18-680.html.
  * Loukas & Vandergheynst (2018) Andreas Loukas and Pierre Vandergheynst.  Spectrally Approximating Large Graphs with Smaller Graphs.  In _Proceedings of the 35th International Conference on Machine Learning_ , pp. 3243–3252. PMLR, 2018.
  * Mavromatis & Karypis (2021) Costas Mavromatis and G. Karypis.  Graph infoclust: Maximizing coarse-grain mutual information in graphs.  In _PAKDD_ , 2021.
  * McInnes et al. (2018) Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger.  Umap: Uniform manifold approximation and projection.  _Journal of Open Source Software_ , 3(29):861, 2018.  doi: 10.21105/joss.00861.  URL https://doi.org/10.21105/joss.00861.
  * Mondal et al. (2021) Anindya Mondal, Shashant R, Jhony H. Giraldo, Thierry Bouwmans, and Ananda S. Chowdhury.  Moving object detection for event-based vision using graph spectral clustering.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV) Workshops_ , pp. 876–884, October 2021.
  * Mrabah et al. (2022) Nairouz Mrabah, Mohamed Bouguessa, Mohamed Fawzi Touati, and Riadh Ksantini.  Rethinking graph auto-encoder models for attributed graph clustering.  _IEEE Transactions on Knowledge and Data Engineering_ , pp. 1–15, 2022.  doi: 10.1109/TKDE.2022.3220948.
  * Newman (2004) M. E. J. Newman.  Fast algorithm for detecting community structure in networks.  _Phys. Rev. E_ , 69:066133, Jun 2004.  doi: 10.1103/PhysRevE.69.066133.  URL https://link.aps.org/doi/10.1103/PhysRevE.69.066133.
  * Newman (2006a) M. E. J. Newman.  Finding community structure in networks using the eigenvectors of matrices.  _Phys. Rev. E_ , 74:036104, Sep 2006a.  doi: 10.1103/PhysRevE.74.036104.  URL https://link.aps.org/doi/10.1103/PhysRevE.74.036104.
  * Newman (2006b) Mark E. J. Newman.  Modularity and community structure in networks.  _Proceedings of the National Academy of Sciences of the United States of America_ , 103 23:8577–82, 2006b.
  * Nowicki & Snijders (2001) Krzysztof Nowicki and Tom A. B Snijders.  Estimation and prediction for stochastic blockstructures.  _Journal of the American Statistical Association_ , 96(455):1077–1087, 2001.  doi: 10.1198/016214501753208735.  URL https://doi.org/10.1198/016214501753208735.
  * Pan et al. (2018) Shirui Pan, Ruiqi Hu, Guodong Long, Jing Jiang, Lina Yao, and Chengqi Zhang.  Adversarially regularized graph autoencoder for graph embedding.  In _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence, IJCAI-18_ , pp. 2609–2615. International Joint Conferences on Artificial Intelligence Organization, 7 2018\.  doi: 10.24963/ijcai.2018/362.  URL https://doi.org/10.24963/ijcai.2018/362.
  * Park et al. (2019) J. Park, M. Lee, H. Chang, K. Lee, and J. Choi.  Symmetric graph convolutional autoencoder for unsupervised graph representation learning.  In _2019 IEEE/CVF International Conference on Computer Vision (ICCV)_ , pp. 6518–6527, Los Alamitos, CA, USA, nov 2019. IEEE Computer Society.  doi: 10.1109/ICCV.2019.00662.  URL https://doi.ieeecomputersociety.org/10.1109/ICCV.2019.00662.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * Perozzi et al. (2014a) Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: online learning of social representations.  In _Proceedings of the 20th ACM SIGKDD international conference on Knowledge discovery and data mining_ , KDD ’14. ACM, August 2014a.  doi: 10.1145/2623330.2623732.  URL http://dx.doi.org/10.1145/2623330.2623732.
  * Perozzi et al. (2014b) Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations.  In _Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’14, pp. 701–710, New York, NY, USA, 2014b. Association for Computing Machinery.  ISBN 9781450329569.  doi: 10.1145/2623330.2623732.  URL https://doi.org/10.1145/2623330.2623732.
  * Rajawat & Kumar (2017) Ketan Rajawat and Sandeep Kumar.  Stochastic multidimensional scaling.  _IEEE Transactions on Signal and Information Processing over Networks_ , 3(2):360–375, 2017.  doi: 10.1109/TSIPN.2017.2668145.
  * Salha-Galvan et al. (2022) Guillaume Salha-Galvan, Johannes F Lutzeyer, George Dasoulas, Romain Hennequin, and Michalis Vazirgiannis.  Modularity-aware graph autoencoders for joint community detection and link prediction.  _Neural Networks_ , 153:474–495, 2022.
  * Shi & Malik (2000) Jianbo Shi and J. Malik.  Normalized cuts and image segmentation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 22(8):888–905, 2000.  doi: 10.1109/34.868688.
  * Tang et al. (2008) Jie Tang, Jing Zhang, Limin Yao, Juanzi Li, Li Zhang, and Zhong Su.  Arnetminer: Extraction and mining of academic social networks.  In _Proceedings of the 14th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’08, pp. 990–998, New York, NY, USA, 2008. Association for Computing Machinery.  ISBN 9781605581934.  doi: 10.1145/1401890.1402008.  URL https://doi.org/10.1145/1401890.1402008.
  * Tsitsulin et al. (2023) Anton Tsitsulin, John Palowitch, Bryan Perozzi, and Emmanuel MÃ¼ller.  Graph clustering with graph neural networks.  _Journal of Machine Learning Research_ , 24(127):1–21, 2023.  URL http://jmlr.org/papers/v24/20-998.html.
  * Veličković et al. (2019) Petar Veličković, William Fedus, William L. Hamilton, Pietro Liò, Yoshua Bengio, and R Devon Hjelm.  Deep Graph Infomax.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=rklz9iAcKQ.
  * Wang et al. (2019) Chun Wang, Shirui Pan, Ruiqi Hu, Guodong Long, Jing Jiang, and Chengqi Zhang.  Attributed graph clustering: A deep attentional embedding approach.  In _Proceedings of the Twenty-Eighth International Joint Conference on Artificial Intelligence, IJCAI-19_ , pp. 3670–3676. International Joint Conferences on Artificial Intelligence Organization, 7 2019\.  doi: 10.24963/ijcai.2019/509.  URL https://doi.org/10.24963/ijcai.2019/509.
  * Wei & Cheng (1989) Yen-Chuen Wei and Chung-Kuan Cheng.  Towards efficient hierarchical designs by ratio cut partitioning.  In _1989 IEEE International Conference on Computer-Aided Design. Digest of Technical Papers_ , pp. 298–301, 1989.  doi: 10.1109/ICCAD.1989.76957.
  * Ying et al. (2018) Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), _Advances in Neural Information Processing Systems_ , volume 31. Curran Associates, Inc., 2018.  URL https://proceedings.neurips.cc/paper_files/paper/2018/file/e77dbaf6759253c7c6d0efc5690369c7-Paper.pdf.
  * Zhang et al. (2019) Xiaotong Zhang, Han Liu, Qimai Li, and Xiao-Ming Wu.  Attributed graph clustering via adaptive graph convolution.  In _Proceedings of the 28th International Joint Conference on Artificial Intelligence_ , IJCAI’19, pp. 4327–4333. AAAI Press, 2019.  ISBN 9780999241141.
  * Zhao et al. (2021) Han Zhao, Xu Yang, Zhenru Wang, Erkun Yang, and Cheng Deng.  Graph debiased contrastive learning with joint representation clustering.  In Zhi-Hua Zhou (ed.), _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21_ , pp. 3434–3440. International Joint Conferences on Artificial Intelligence Organization, 8 2021.  doi: 10.24963/ijcai.2021/473.  URL https://doi.org/10.24963/ijcai.2021/473.  Main Track.
  * Zhao et al. (2012) Yunpeng Zhao, Elizaveta Levina, and Ji Zhu.  Consistency of community detection in networks under degree-corrected stochastic block models.  _The Annals of Statistics_ , 40(4):2266 – 2292, 2012.  doi: 10.1214/12-AOS1036.  URL https://doi.org/10.1214/12-AOS1036.

Appendix

[sections] [sections]l1

##  Appendix A Convexity of terms in the optimization objective 4

When $X_{C}$ is kept constant, $\mathcal{L}_{MAGC}$ gets reduced to:

$\displaystyle\underset{C}{\min}\;f(C)=\ \textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})+$
$\displaystyle\frac{\alpha}{2}\left\|CX_{C}-X\right\|_{F}^{2}-\frac{\beta}{2e}\textrm{tr}(C^{T}BC)$
(24) $\displaystyle-\gamma\log\det(C^{T}\Theta
C+J)+\frac{\lambda}{2}\left\|C^{T}\right\|_{1,2}^{2}$
$\displaystyle\text{subject}\;\;\text{to}\
\mathcal{C}\in\mathcal{S}_{c}\ref{coarsened_features_equation}\;\text{where,
}J=\frac{1}{k}\textbf{1}_{k\times k}$

The term $\textrm{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$ is convex function in $C$.
This result can be derived easily using Cholesky Decomposition on the positive
semi-definite matrix $\Theta$ (i.e. $\Theta=L^{T}L$):

$\displaystyle\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})=\textrm{tr}(X_{C}^{T}C^{T}L^{T}LCX_{C})=\textrm{tr}(({LCX_{C}})^{T}LCX_{C})=\left\|LCX_{C}\right\|_{F}^{2}$
(26)

Frobenius norm is a convex function, and the simplified expression is linear
in C. Hence we can deduce that the tr(.) term is convex in C. The terms
$\left\|CX_{C}-X\right\|_{F}^{2}$ and $\left\|C^{T}\right\|_{1,2}^{2}$ are
convex because Frobenius norm and $l_{1,2}$ norm are convex in C.

For proving the convexity of $-\log\det(C^{T}\Theta C+J)$ we restrict function
to a line. We define a function $g$:

$\displaystyle g(t)=f(z+tu)\text{where, }t\in dom(g),z\in
dom(f),u\in\mathbb{R}^{n}.$ (27)

A function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is convex if
$g:\mathbb{R}\rightarrow\mathbb{R}$ is convex.

The graph Laplacian matrix of the coarsened graph ($\Theta_{c}$) is symmetric
and positive semi-definite having a rank of k-1. To convert $\Theta_{c}$ to
positive definite matrix, we add a rank 1 matrix
$J=\frac{1}{k}\textbf{1}_{k\times k}$. ($\Theta_{c}+J=L^{T}L$)

$\displaystyle f(L)=-\log\det(C^{T}\Theta C+J)=-\log\det(L^{T}L)$ (28)

Now substituting L = Z + tU in the above equation.

$\displaystyle g(t)$ $\displaystyle=-\log\det((Z+tU)^{T}(Z+tU))$ (29)
$\displaystyle=-\log\det(Z^{T}Z+t(Z^{T}U+U^{T}Z)t^{2}U^{T}U)$ (30)
$\displaystyle=-\log\det(Z^{T}(I+t(UZ^{-1}+{(UZ^{-1})}^{T})+t^{2}{(Z^{-1})}^{T}U^{T}UZ^{-1})Z)$
(31) $\displaystyle\text{substituting }\ P$ $\displaystyle=VZ^{-1}$ (32)
$\displaystyle=-(\log\det(Z^{T}Z)+\log\det(I+t(P+P^{T})+t^{2}P^{T}P))$ (33)
$\displaystyle\text{Eigenvalue decomposition of}P$ $\displaystyle=Q\Lambda
Q^{T}\text{and}QQ^{T}=I$ (34)
$\displaystyle=-(\log\det(Z^{T}Z)+\log\det(QQ^{T}+2tQ\Lambda
Q^{T}+t^{2}Q\Lambda^{2}Q^{T}))$ (35)
$\displaystyle=-(\log\det(Z^{T}Z)+\log\det(Q(I+2t\Lambda+t^{2}\Lambda^{2})Q^{T}))$
(36)
$\displaystyle=-\log\det(Z^{T}Z)-\sum_{i=1}^{n}\log(1+2t\lambda_{i}+t^{2}\lambda^{2})$
(37)

Finding double derivative of $g(t)$:

$\displaystyle
g"(t)=\sum_{i=1}^{n}\frac{2\lambda_{i}^{2}(1+t\lambda_{i})^{2}}{(1+2t\lambda_{i}+t^{2}\lambda_{i}^{2})^{2}}$
(38)

Since $g"(t)\geq 0\forall\ t\in\mathbb{R}$, $g(t)$ is a convex function in
$t$. This implies $f(L)$ is convex in $L$. We know that, $C^{T}\Theta
C+J=L^{T}L$ so,

$\displaystyle L=\Theta^{\frac{1}{2}}C+\frac{1}{\sqrt{kp}}\textbf{1}_{p\times
k}$ (39)

Since $L$ is linear in $C$ and $f(L)$ is convex in $L$, $-\log\det(C^{T}\Theta
C+J)$ is convex in C.

##  Appendix B Optimal Solution of Optimization Objective in Equation 4

We first show that the function $f(C)$ is $L-Lipschitz$ continuous gradient
function with $L=\max(L_{1},L_{2},L_{3},L_{4},L_{5})$, where
$L_{1},L_{2},L_{3},L_{4},andL_{5}$ are the Lipschitz constants of
$\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C}),\frac{\alpha}{2}\left\|CX_{C}-X\right\|_{F}^{2},-\frac{\beta}{2e}\textrm{tr}(C^{T}BC),-\gamma\log\det(C^{T}\Theta
C+J),and\frac{\lambda}{2}\left\|C^{T}\right\|_{1,2}^{2}$.

For the $\textrm{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$ term, we apply triangle
inequality and employ the property of the norm of the trace operator:
$||tr||=\underset{M\neq 0}{sup}\frac{|tr(M)|}{||M||_{F}}$.

$\displaystyle|tr(X_{C}^{T}C_{1}^{T}\Theta
C_{1}X_{C})-tr(X_{C}^{T}C_{2}^{T}\Theta C_{2}X_{C})|$ (40)
$\displaystyle=|tr(X_{C}^{T}C_{1}^{T}\Theta
C_{1}X_{C})-tr(X_{C}^{T}C_{2}^{T}\Theta
C_{1}X_{C})+tr(X_{C}^{T}C_{2}^{T}\Theta
C_{1}X_{C})-tr(X_{C}^{T}C_{2}^{T}\Theta C_{2}X_{C})|$ (41)
$\displaystyle\leq|tr(X_{C}^{T}C_{1}^{T}\Theta
C_{1}X_{C})-tr(X_{C}^{T}C_{2}^{T}\Theta
C_{1}X_{C})|+|tr(X_{C}^{T}C_{2}^{T}\Theta
C_{1}X_{C})-tr(X_{C}^{T}C_{2}^{T}\Theta C_{2}X_{C})|$ (42)
$\displaystyle\leq||tr||||X_{C}^{T}(C_{1}-C_{2})^{T}\Theta
C_{1}X_{C}||_{F}+||tr||||X_{C}^{T}C_{2}^{T}\Theta(C_{1}-C_{2})X_{C}||_{F}$
(43)
$\displaystyle\leq||tr||||X_{C}||_{F}||\Theta||||C_{1}-C_{2}||_{F}(||C_{1}||_{F}+||C_{2}||_{F})\
\ \text{(Frobenius Norm Property)}$ (44) $\displaystyle\leq
2\sqrt{p}||tr||||X_{C}||_{F}||\Theta||||C_{1}-C_{2}||_{F}\ \
(||C_{1}||_{F}=||C_{2}||_{F}=\sqrt{p})$ (45) $\displaystyle\leq
L_{1}||C_{1}-C_{2}||_{F}$ (46)

The second term is $\frac{\alpha}{2}\left\|CX_{C}-X\right\|_{F}^{2}$ can be
written as:

$\displaystyle\frac{\alpha}{2}tr((CX_{C}-X)^{T}(CX_{C}-X))$ (47)
$\displaystyle=\frac{\alpha}{2}tr(X_{C}^{T}C^{T}CX_{C}-X^{T}CX_{C}+X^{T}X-X_{C}^{T}C^{T}X)$
(48)
$\displaystyle=\frac{\alpha}{2}(tr(X_{C}^{T}C^{T}CX_{C})-tr(X^{T}CX_{C})+tr(X^{T}X)-tr(X_{C}^{T}C^{T}X))$
(49)

All the terms except $tr(X^{T}X)$ (constant with respect to C) in obtained in
the expression will follow similar proofs to $\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})$.

Next we consider the modularity term:

$\displaystyle|tr(C_{1}^{T}BC_{1})-tr(C_{2}^{T}BC_{2})|$ (50)
$\displaystyle=|tr(C_{1}^{T}BC_{1})-tr(C_{2}^{T}BC_{1})+tr(C_{2}^{T}BC_{1})-tr(C_{2}^{T}BC_{2})|$
(51)
$\displaystyle\leq|tr(C_{1}^{T}BC_{1})-tr(C_{2}^{T}BC_{1})|+|tr(C_{2}^{T}BC_{1})-tr(C_{2}^{T}BC_{2})|$
(52)
$\displaystyle\leq||tr||||(C_{1}-C_{2})^{T}BC_{1}||_{F}+||tr||||(C_{1}-C_{2})^{T}BC_{2}||_{F}$
(53)
$\displaystyle\leq||tr||||B||||C_{1}-C_{2}||_{F}(||C_{1}||_{F}+||C_{2}||_{F})\
\ \text{(Frobenius Norm Property)}$ (54) $\displaystyle\leq
L_{3}||C_{1}-C_{2}||_{F}$ (55)

The Lipschitz constant for $-\gamma\log\det(C^{T}\Theta C+J)$ is linked to the
smallest non-zero eigenvalue of the coarsened Laplacian matrix ($\Theta_{c}$)
and is bounded by $\frac{\delta}{(k-1)^{2}}$ (Rajawat & Kumar, 2017), where
$\delta$ is the minimum non-zero weight of $G_{c}$.

$tr(\textbf{1}^{T}C^{T}C\textbf{1})$ can be proved to be $L_{5}-Lipschitz$
like the modularity and Dirichlet energy (smoothness) terms. This concludes
the proof.

The majorized problem for L-Lipschitz and differentiable functions can now be
applied. The Lagrangian of the majorized problem, 8 is:

$\displaystyle\mathcal{L}(C,X_{C},\mu)=\frac{1}{2}C^{T}C-C^{T}A-\mu_{1}^{T}C+\mu_{2}^{T}\Big{[}\left\|C_{1}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{i}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{p}^{T}\right\|_{2}^{2}-1\Big{]}^{T}$
(56)

where $\mu=\mu_{1}||\mu_{2}$ are the dual variables and
$A=\Big{(}C-\frac{1}{L}\nabla f(C)\Big{)}^{+}$

The corresponding KKT conditions (w.r.t $C$) are:

$\displaystyle
C-A-\mu_{1}+2[{\mu_{2}}_{o}C_{0}^{T},\cdots,{\mu_{2}}_{i}C_{i}^{T},\cdots,{\mu_{2}}_{p}C_{p}^{T}]$
$\displaystyle=0$ (57)
$\displaystyle\mu_{2}^{T}\Big{[}\left\|C_{1}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{i}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{p}^{T}\right\|_{2}^{2}-1,\Big{]}^{T}$
$\displaystyle=0$ (58) $\displaystyle\mu_{1}^{T}C$ $\displaystyle=0$ (59)
$\displaystyle\mu_{1}$ $\displaystyle\geq 0$ (60) $\displaystyle\mu_{2}$
$\displaystyle\geq 0C$ $\displaystyle\geq 0$ (61)
$\displaystyle\left\|[C^{T}]_{i}\right\|_{2}^{2}$ $\displaystyle\leq 1\
\forall i$ (62)

The optimal solution to these KKT conditions is:

$\displaystyle C=\frac{(A)^{+}}{\sum_{i}\left\|[A^{T}]_{i}\right\|}_{2}$ (63)

##  Appendix C Proof of Theorem 1 (Convergence)

In this section, we prove that the sequence $\\{C^{t+1},X_{C}^{t+1}\\}$
generated by Algorithm 1 converges to the set of Karush–Kuhn–Tucker (KKT)
optimality points for Problem 4.

The Lagrangian of Problem 4 comes out to be:

$\displaystyle\mathcal{L}(C,X_{C},\mathbf{\mu})$
$\displaystyle=\textrm{tr}(X_{C}^{T}C^{T}\Theta
CX_{C})+\frac{\alpha}{2}\left\|CX_{C}-X\right\|_{F}^{2}-\frac{\beta}{2e}\textrm{tr}(C^{T}BC)$
(64) $\displaystyle\ \ -\gamma\log\det(C^{T}\Theta
C+J)+\frac{\lambda}{2}\left\|C^{T}\right\|_{1,2}^{2}-\mathbf{\mu}_{1}^{T}C+\sum_{i}{\mathbf{\mu}_{2}}_{i}\Big{[}\left\|C^{T}_{i}\right\|_{2}^{2}-1\Big{]}$
(65)

where $\mu=\mu_{1}||\mu_{2}$ are the dual variables.

w.r.t. $C$, the KKT conditions are

$\displaystyle 2\Theta
CX_{C}X_{C}^{T}+\alpha(CX_{C}-X)X_{C}^{T}-\frac{\beta}{e}BC-2\gamma\Theta
C(C^{T}\Theta C+J)^{-1}$ (66) $\displaystyle+\lambda C\mathbf{1}_{k\times
k}-\mu_{1}+2[{\mu_{2}}_{o}C_{0}^{T},\cdots,{\mu_{2}}_{i}C_{i}^{T},\cdots,{\mu_{2}}_{p}C_{p}^{T}]$
$\displaystyle=0$
$\displaystyle\mu_{2}^{T}\Big{[}\left\|C_{1}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{i}^{T}\right\|_{2}^{2}-1,\cdots,\left\|C_{p}^{T}\right\|_{2}^{2}-1\Big{]}^{T}$
$\displaystyle=0$ (67) $\displaystyle\mu_{1}^{T}C$ $\displaystyle=0$ (68)
$\displaystyle\mu_{1}$ $\displaystyle\geq 0$ (69) $\displaystyle\mu_{2}$
$\displaystyle\geq 0$ (70) $\displaystyle C$ $\displaystyle\geq 0$ (71)
$\displaystyle\left\|[C^{T}]_{i}\right\|_{2}^{2}$ $\displaystyle\leq 1\
\forall i$ (72)

Now, $C^{\infty}\equiv\underset{t\rightarrow\infty}{\lim}{C^{t}}$ is found
from Equation 9 as:

$\displaystyle C^{\infty}$
$\displaystyle=C^{\infty}+\frac{1}{L}\Bigg{(}2\Theta
C^{\infty}X_{C}^{\infty}X_{C}^{\infty}+\alpha(C^{\infty}X_{C}-X)X_{C}^{\infty}-\frac{\beta}{e}BC^{\infty}$
(73) $\displaystyle\ \ \ \ -2\gamma\Theta C^{\infty}({C^{\infty}}^{T}\Theta
C^{\infty}+J)^{-1}+\lambda C^{\infty}\mathbf{1}_{k\times k}\Bigg{)}$
$\displaystyle 0$ $\displaystyle=2\Theta
C^{\infty}X_{C}^{\infty}X_{C}^{\infty}+\alpha(C^{\infty}X_{C}-X)X_{C}^{\infty}-\frac{\beta}{e}BC^{\infty}$
(74) $\displaystyle\ \ \ \ -2\gamma\Theta C^{\infty}({C^{\infty}}^{T}\Theta
C^{\infty}+J)^{-1}+\lambda C^{\infty}\mathbf{1}_{k\times k}$

So, for $\mu=0$, $C^{\infty}$ satisfies the KKT conditions.

w.r.t. $X_{C}$, the KKT conditions are:

$\displaystyle 2C^{T}\Theta CX_{C}+\alpha C^{T}(CX_{C}-X)=0$ (75)

So, $X^{\infty}\equiv\underset{t\rightarrow\infty}{\lim}{X^{t}}$ found from
Equation 13 will satisfy this as that equation is just a rearrangement of the
KKT condition.

##  Appendix D Proof (continued) of Theorem 2 (conditions for consistency)

As previously defined in Theorem 2 $O=C^{T}AC$ and

$\frac{1}{\mu_{p}}\mathbb{E}[O|C,t]=H(S)\qquad$

We need to find "population version" of the loss function in terms of $H(S)$.

$\displaystyle\mathbb{E}[O|C,t]$
$\displaystyle=\mathbb{E}[C^{T}AC|C,t]=\mu_{p}H(S)$ (76)
$\displaystyle\mathbb{E}[C^{T}DC|C,t]$
$\displaystyle=\mathbb{E}[C^{T}\text{diag}(\sum_{j=1}^{k}A_{ij})C|C,t]$
$\displaystyle=\mathbb{E}[\text{diag}(\sum_{j=1}^{k}O_{ij})|C,t]$ (77)
$\displaystyle=\mu_{p}\text{diag}(\sum_{j=1}^{k}H_{ij})$ (78)

So, for $\text{tr}(X_{C}^{T}C^{T}\Theta CX_{C})$

$\displaystyle\mathbb{E}[\text{tr}(X_{C}^{T}C^{T}\Theta CX_{C})|C,t]$
$\displaystyle=\text{tr}(\mathbb{E}[X_{C}^{T}C^{T}\Theta CX_{C}|C,t])$
$\displaystyle=\text{tr}(X_{C}^{T}\mathbb{E}[C^{T}\Theta C|C,t]X_{C})$
$\displaystyle=\text{tr}(X_{C}^{T}\mathbb{E}[C^{T}DC-C^{T}AC|C,t]X_{C})$
$\displaystyle=\text{tr}(X_{C}^{T}[\mu_{p}\text{diag}(\sum_{j=1}^{k}H_{ij})$
$\displaystyle-\mu_{p}H(S)]X_{C})$
$\displaystyle=\mu_{p}\text{tr}(X_{C}^{T}[\text{diag}(\sum_{j=1}^{k}H_{ij})$
$\displaystyle-H(S)]X_{C})$ (79)

Next, we have $\frac{1}{2e}\text{tr}(C^{T}BC)$, which has already been solved
in the paper Zhao et al. (2012) in their Appendix (Page 25 of the full
document).

$\frac{1}{2e}\mathbb{E}[\text{tr}(C^{T}BC)]=\sum_{k}\bigg{(}\frac{H_{kk}}{\tilde{P_{0}}}-\bigg{(}\frac{H_{k}}{\tilde{P_{0}}}\bigg{)}^{2}\bigg{)}$

For $\log\det(C^{T}\Theta C+J)$, where $J=\frac{1}{k}\mathbb{1}_{k\times k}$,

We can write $\log(\det(C^{T}\Theta C+J))=\text{tr}(\log(C^{T}\Theta C+J))$
since, $\det(A)=e^{\text{tr}(\log(A))}$.

$\displaystyle Z=C^{T}\Theta C+J$ $\displaystyle=V\Lambda V^{-1}$ (80) So,
(81) $\displaystyle\text{tr}(\log(V\Lambda V^{-1}))$
$\displaystyle=\text{tr}(V\log(\Lambda)V^{-1})$ (82)
$\displaystyle\log\Lambda$
$\displaystyle=\log(bI)+\log\bigg{(}I+\frac{\Lambda}{b}-I\bigg{)}$ (83)

Using the first-order Taylor expansion of $\log(I+X)=X$, we need to choose $b$
such that

$l=\bigg{\|}\frac{\Lambda}{b}-I\bigg{\|}_{F}<1$ (84)

And for the expansion to be a good approximation, we need $l\to 0$. We will
enforce this later in equation 94.

$\displaystyle\text{tr}(V\bigg{(}\log(bI)+\log\bigg{(}I+\frac{\Lambda}{b}-I\bigg{)}\bigg{)}V^{-1})$
(85) $\displaystyle=$ $\displaystyle\
\text{tr}(V\bigg{(}\log(bI)+\bigg{(}I+\frac{\Lambda}{b}-I\bigg{)}\bigg{)}V^{-1})$
(86) $\displaystyle=$ $\displaystyle\
\text{tr}(V\log(bI)V^{-1})+\frac{1}{b}\text{tr}(V\Lambda V^{-1})-\text{tr}(I)$
(87) $\displaystyle=$ $\displaystyle\
\text{tr}(\log(bI))+\frac{1}{b}\text{tr}(Z)-\text{tr}(I)$ (88)
$\displaystyle=$ $\displaystyle\ k\log(b)+\frac{1}{b}\text{tr}(Z)-k$ (89)

Finding the expectation of from equation 83,

$\displaystyle=k\log(b)+\mathbb{E}\bigg{[}\frac{1}{b}\text{tr}(Z)\bigg{|}\
c,t\bigg{]}-k$ (90)
$\displaystyle=k\log(b)+\frac{1}{b}\text{tr}(\mathbb{E}[Z|\ c,t])-k$ (91)
$\displaystyle=k\log(b)+\frac{1}{b}\text{tr}(\mathbb{E}[C^{T}\Theta C+J|\
c,t])-k$ (92) $\displaystyle\text{Using $\Theta=D-A$,
equation~{}\ref{eqn:supp_expect_CTAC} and
equation~{}\ref{eqn:supp_expect_CTDC}},$
$\displaystyle=k\log(b)+\frac{1}{b}\text{tr}\bigg{(}\mu_{p}\text{diag}(\sum_{j=1}^{k}H_{ij})-\mu_{p}H(S)\bigg{)}-k$
(93)

which is a linear function in $H(S)$.

For the approximation to be good, we can now simplify $l$ as defined in
equation 84:

$\displaystyle l$ $\displaystyle=\bigg{\|}\frac{\Lambda}{b}-I\bigg{\|}_{F}$
(94) $\displaystyle l$
$\displaystyle=\begin{Vmatrix}\frac{1}{b}\Lambda_{11}-1&\frac{1}{b}\Lambda_{12}&\cdots&\frac{1}{b}\Lambda_{1k}\\\
\frac{1}{b}\Lambda_{21}&\frac{1}{b}\Lambda_{22}-1&\cdots&\frac{1}{b}\Lambda_{2k}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{1}{b}\Lambda_{k1}&\frac{1}{b}\Lambda_{k2}&\cdots&\frac{1}{b}\Lambda_{kk}-1\\\
\end{Vmatrix}_{F}$ (95)

Writing out this norm, we get a quadratic expression in $\frac{1}{b}$:

$\displaystyle l$
$\displaystyle=(\sum_{i=1}^{k}\sum_{j=1}^{k}\lambda_{ij}^{2})\frac{1}{b^{2}}-2(\sum_{u=1}^{k}\lambda_{ii})\frac{1}{b}+k$
(96) $\displaystyle\text{or concisely},l$
$\displaystyle=\|\Lambda\|_{F}^{2}\frac{1}{b^{2}}-2\text{tr}(\Lambda)\frac{1}{b}+k$
(97) $\displaystyle\text{Since, }\Lambda$ is a diagonal matrix, (98)
$\displaystyle\text{tr}(\Lambda)$
$\displaystyle=\sum_{i}\lambda_{i}=\text{tr}(Z)$ (99) Also, since
$\displaystyle\Lambda^{2}\text{ is the eigenvalue matrix for }Z^{2},$ (100)
$\displaystyle\|\Lambda\|_{F}^{2}$
$\displaystyle=\sum_{i}\lambda_{i}^{2}=\text{tr}(Z^{2})$ (101) $\displaystyle
l$ $\displaystyle=\text{tr}(Z^{2})\frac{1}{b^{2}}-2\text{tr}(Z)\frac{1}{b}+k$
(102)

Since, $l<1\implies l-1<0\implies l-1=0$ has 2 real roots. Using simple
quadratic analysis (in $\frac{1}{b}$), the discriminant $\Delta$ should be
positive.

$\displaystyle\Delta=4\text{tr}(Z)^{2}\ -4(k-1)\text{tr}(Z^{2})$
$\displaystyle>0$ (103)
$\displaystyle\frac{\text{tr}(Z)^{2}}{\text{tr}(Z^{2})}$ $\displaystyle>k-1$
(104)

The minimum $\displaystyle\text{ value of }l\coloneqq l_{min}\text{ occurs at
}$ (105) $\displaystyle b$
$\displaystyle=\frac{\text{tr}(Z^{2})}{\text{tr}(Z)}=\text{tr}(Z)\cdot\frac{\text{tr}(Z^{2})}{\text{tr}(Z)^{2}}$
(106) $\displaystyle l_{min}$ $\displaystyle<1\text{ will exist when
}\frac{\text{tr}(Z)^{2}}{\text{tr}(Z^{2})}>k-1$ (107) which holds for
$\displaystyle b<\frac{\text{tr}(Z)}{k-1}$ (108) $\displaystyle\text{So we can
always choose }b$ $\displaystyle<\frac{2k-1}{k-1}$ (109) $\displaystyle\text{
since }\min\text{tr}(Z)$ $\displaystyle=2k-1[\text{proved in
\ref{eqn:supp_min_trZ}}]$ (110)

$\displaystyle\min\text{tr}(C^{T}\Theta C+J)$
$\displaystyle=\min\text{tr}(C^{T}\Theta C)+\text{tr}(J)$ (111)
$\displaystyle=\min\text{tr}(C^{T}DC)-\text{tr}(C^{T}AC)+1$ (112)
$\displaystyle=2e-2(e-(k-1))+1$ (113) $\displaystyle=2k-1$ (114)

Additionally, define $\tilde{\pi}_{a}=\sum_{u}x_{u}\Pi_{au}$ with
$\sum_{a}\tilde{\pi}_{a}=1$, since $\mathbb{E}[t_{i}]=1$.

###  D.1 Required Condition a): Lipschitz Continuity

Condition 1: We need to show that
$|F(S_{1})-F(S_{2})|\leq\alpha\|S_{1}-S_{2}\|$ (Lipschitz)

$\displaystyle|F(S_{1})-F(S_{2})|$
$\displaystyle\leq|f_{1}(S_{1})-f_{1}(S_{2})|+|f_{2}(S_{1})-f_{2}(S_{2})|+|f_{3}(S_{1})-f_{3}(S_{2})|$
(115)

Let’s first find $\|H(S_{1})-H(S_{2})\|_{F}$

$\displaystyle H(S)$ $\displaystyle=(Sx)P(Sx)^{T}$ (116)
$\displaystyle\|H(S_{1})-H(S_{2})\|_{F}$
$\displaystyle=\|(S_{1}x)P(S_{1}x)^{T}-(S_{2}x)P(S_{2}x)^{T}\|_{F}$ (117)
$\displaystyle=\|(S_{1}x)P((S_{1}-S_{2})x)^{T}+((S_{1}-S_{2})x)P(S_{2}x)^{T}\|_{F}$
(118)

Next, we see
$\|\text{diag}(\sum_{j=1}^{k}H(S_{1})_{ij})-\text{diag}(\sum_{j=1}^{k}H(S_{2})_{ij})\|_{F}=\|\text{diag}(\sum_{j=1}^{k}(H(S_{1})_{ij}-H(S_{2})_{ij}))\|_{F}$

For the first term, $|f_{1}(H(S_{1}))-f_{1}(H(S_{2}))|$, define
$H(S)_{i}=\sum_{j=1}^{k}H(S)_{ij}$

$\displaystyle=|\mu_{p}\text{tr}(X_{C}^{T}[H(S_{1})-\text{diag}(\sum_{j=1}^{k}H(S_{1})_{ij})]X_{C})-\mu_{p}\text{tr}(X_{C}^{T}[H(S_{2})-\text{diag}(\sum_{j=1}^{k}H(S_{2})_{ij})]X_{C})|$
(119)
$\displaystyle=|\mu_{p}\bigg{(}\text{tr}(X_{C}^{T}[H(S_{1})-H(S_{2})]X_{C})-\text{tr}(X_{C}^{T}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])X_{C})\bigg{)}|$
(120)
$\displaystyle\leq|\mu_{p}|\bigg{(}\bigg{|}\text{tr}(X_{C}^{T}[H(S_{1})-H(S_{2})]X_{C})\bigg{|}-\bigg{|}\text{tr}(X_{C}^{T}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])X_{C})\bigg{|}\bigg{)}$
(121)
$\displaystyle\leq|\mu_{p}|\bigg{(}\|\text{tr}\|\bigg{\|}X_{C}^{T}[H(S_{1})-H(S_{2})]X_{C}\bigg{\|}_{F}-\|\text{tr}\|\bigg{\|}X_{C}^{T}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])X_{C}\bigg{\|}_{F}\bigg{)}$
(122)
$\displaystyle\leq|\mu_{p}|\bigg{(}\|\text{tr}\|\|X_{C}\|^{2}\bigg{\|}H(S_{1})-H(S_{2})\bigg{\|}_{F}-\|\text{tr}\|\|X_{C}\|^{2}\bigg{\|}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])\bigg{\|}_{F}\bigg{)}$
(123)

Taking the first sub-term,

$\displaystyle|\mu_{p}|\ \|\text{tr}\|\ \|X_{C}\|^{2}\
\bigg{\|}H(S_{1})-H(S_{2})\bigg{\|}_{F}$ (124) $\displaystyle=\|\text{tr}\|\
\|X_{C}\|^{2}\
\bigg{\|}(S_{1}x)P((S_{1}-S_{2})x)^{T}+((S_{1}-S_{2})x)P(S_{2}x)^{T}\bigg{\|}_{F}$
(125) $\displaystyle\leq\|\text{tr}\|\ \|X_{C}\|^{2}\ \|P\|_{F}\
(\|S_{1}x\|_{F}+\|S_{2}x\|_{F})\ \|(S_{1}-S_{2})x\|_{F}$ (126)
$\displaystyle\leq\|\text{tr}\|\ \|X_{C}\|^{2}\ \|P\|_{F}\
(\|S_{1}x\|_{F}+\|S_{2}x\|_{F})\ \|t\|_{F}\ \|(S_{1}-S_{2})\|_{F}$ (127)
$\displaystyle=\alpha_{1}\|(S_{1}-S_{2})\|_{F}$ (128)

For the second sub-term, define
$\tilde{S}_{ka}=\sum_{u}x_{u}S_{kau}=(Sx)_{ka}$

$\displaystyle H(S)_{i}=\sum_{j=1}^{k}H(S)_{ij}$
$\displaystyle=\sum_{as}\tilde{\pi}_{s}P_{as}\tilde{S}_{ia}=\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}S_{iau}$
(129) $\displaystyle|H(S_{1})_{i}-H(S_{2})_{i}|$
$\displaystyle=|\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}{S_{1}}_{iau}-\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}{S_{2}}_{iau}|$
(130)
$\displaystyle=|\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}({S_{1}-S_{2}})_{iau}|$
(131)
$\displaystyle\leq\sum_{asu}|\tilde{\pi}_{s}P_{as}x_{u}({S_{1}-S_{2}})_{iau}|$
(132)
$\displaystyle\leq\sum_{asu}|\tilde{\pi}_{s}P_{as}x_{u}|\cdot|({S_{1}-S_{2}})_{iau}|$
(133)
$\displaystyle\leq\sum_{asu}|\tilde{\pi}_{s}P_{as}x_{u}|\cdot\sum_{au}|({S_{1}-S_{2}})_{iau}|$
(134) $\displaystyle=\alpha_{1}^{\prime}\sum_{au}|({S_{1}-S_{2}})_{iau}|$
(135) $\displaystyle\text{So,
}\|\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])\|_{F}$
$\displaystyle=\sqrt{\sum_{i=1}^{k}\bigg{(}\sum_{j=1}^{k}H(S_{1})_{ij}-H(S_{2})_{ij}\bigg{)}^{2}}$
(136)
$\displaystyle\leq\sqrt{\sum_{i=1}^{k}\bigg{(}\alpha_{1}^{\prime}\sum_{au}|({S_{1}-S_{2}})_{iau}|\bigg{)}^{2}}$
(137) $\displaystyle=\alpha_{1}^{\prime}\|S_{1}-S_{2}\|_{1,1,2}$ (138)

For the second term, $|f_{2}(H(S_{1}))-f_{2}(H(S_{2}))|$

$\displaystyle=\frac{\mu_{p}}{b}\bigg{|}\text{tr}(H(S_{1})-\text{diag}([H(S_{1})_{i}]))$
$\displaystyle-\text{tr}(H(S_{2})-\text{diag}([H(S_{2})_{i}]))\bigg{|}$ (139)
$\displaystyle=\frac{\mu_{p}}{b}\bigg{|}\text{tr}(H(S_{1})-H(S_{2}))$
$\displaystyle-\text{tr}\bigg{(}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])\bigg{)}\bigg{|}$
(140)
$\displaystyle\leq\frac{\mu_{p}}{b}\bigg{(}\bigg{|}\text{tr}(H(S_{1})-H(S_{2}))\bigg{|}$
$\displaystyle+\bigg{|}\text{tr}\bigg{(}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])\bigg{)}\bigg{|}\bigg{)}$
(141) $\displaystyle\leq\frac{\mu_{p}}{b}\bigg{(}\|\text{tr}\|\
\bigg{\|}H(S_{1})-H(S_{2})\bigg{\|}_{F}$
$\displaystyle+\|\text{tr}\|\bigg{\|}\text{diag}([H(S_{1})_{i}-H(S_{2})_{i}])\bigg{\|}_{F}\bigg{|}\bigg{)}$
(142) As shown above, (143) $\displaystyle\leq\alpha_{2}\|(S_{1}-S_{2})\|_{F}$
$\displaystyle+\alpha_{2}^{\prime}\|(S_{1}-S_{2})\|_{1,1,2}$ (144)

For the third term, it has already been proven in Zhao et al. (2012), but we
also prove it here:

$|f_{3}(H(S_{1}))-f_{3}(H(S_{2}))|$

$\displaystyle=\bigg{|}\frac{\text{tr}(H(S_{1}))}{\tilde{P_{0}}}-\frac{\sum_{i=1}^{k}(\sum_{j=1}^{k}H(S_{1})_{ij})^{2}}{\tilde{P_{0}}^{2}}-\frac{\text{tr}(H(S_{2}))}{\tilde{P_{0}}}+\frac{\sum_{i=1}^{k}(\sum_{j=1}^{k}H(S_{2})_{ij})^{2}}{\tilde{P_{0}}^{2}}\bigg{|}$
(145)
$\displaystyle=\bigg{|}\frac{\text{tr}(H(S_{1})-H(S_{2}))}{\tilde{P_{0}}}-\frac{\sum_{i=1}^{k}\big{[}H(S_{1})_{i}^{2}-H(S_{2})_{i}^{2}\big{]}}{\tilde{P_{0}}^{2}}\bigg{|}$
(146)
$\displaystyle\leq\frac{1}{\tilde{P_{0}}}\bigg{|}\text{tr}(H(S_{1})-H(S_{2}))\bigg{|}+\frac{1}{\tilde{P_{0}}^{2}}\bigg{|}\sum_{i=1}^{k}\big{[}H(S_{1})_{i}^{2}-H(S_{2})_{i}^{2}\big{]}\bigg{|}$
(147) As shown above, (148)
$\displaystyle\leq\alpha_{3}\|(S_{1}-S_{2})\|_{F}+\frac{1}{\tilde{P_{0}}^{2}}\bigg{|}\sum_{i=1}^{k}\big{[}H(S_{1})_{i}-H(S_{2})_{i}\big{]}\cdot\big{[}H(S_{1})_{i}+H(S_{2})_{i}\big{]}\bigg{|}$
(149)
$\displaystyle\leq\alpha_{3}\|(S_{1}-S_{2})\|_{F}+\frac{1}{\tilde{P_{0}}^{2}}\sum_{i=1}^{k}\big{|}H(S_{1})_{i}+H(S_{2})_{i}\big{|}\cdot\sum_{i=1}^{k}\big{|}H(S_{1})_{i}-H(S_{2})_{i}\big{|}$
(150)
$\displaystyle=\alpha_{3}\|(S_{1}-S_{2})\|_{F}+\alpha_{3}^{\prime}\sum_{iau}\big{|}S_{1}-S_{2}\big{|}_{iau}$
(151)
$\displaystyle=\alpha_{3}\|(S_{1}-S_{2})\|_{F}+\alpha_{3}^{\prime}\|S_{1}-S_{2}\|_{1,1,1}$
(152)

###  D.2 Required Condition b): Continuity of directional second derivative

Condition 2: $W=H(\mathbb{D})$

$\displaystyle\frac{\partial^{2}}{\partial\varepsilon^{2}}F(M_{0}+\varepsilon(M_{1}-M_{0}),\mathbf{t_{0}}+\varepsilon(\mathbf{t_{1}}-\mathbf{t_{0}}))\bigg{|}_{\varepsilon=0^{+}}$
(153)
$\displaystyle=\frac{\partial^{2}}{\partial\varepsilon^{2}}f_{1}(M_{0}+\varepsilon(M_{1}-M_{0}))+\frac{\partial^{2}}{\partial\varepsilon^{2}}f_{2}(M_{0}+\varepsilon(M_{1}-M_{0}))+\frac{\partial^{2}}{\partial\varepsilon^{2}}f_{3}(M_{0}+\varepsilon(M_{1}-M_{0}))\bigg{|}_{\varepsilon=0^{+}}$
(154)

Finding the directional derivative for $f_{1}$,

$\displaystyle f_{1}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=\mu_{p}\text{tr}(X_{C}^{T}[M_{0}+\varepsilon(M_{1}-M_{0})-\text{diag}(\sum_{j=1}^{k}\big{(}M_{0}+\varepsilon(M_{1}-M_{0})\big{)}_{ij})]X_{C})$
(155)
$\displaystyle=\mu_{p}\bigg{(}\text{tr}(X_{C}^{T}M_{0}X_{C})+\varepsilon\
\text{tr}(X_{C}^{T}(M_{1}-M_{0})X_{C})$ (156)
$\displaystyle\quad-\text{tr}(X_{C}^{T}\text{diag}([\sum_{j=1}^{k}(M_{0})_{ij}])X_{C})-\varepsilon\
\text{tr}(X_{C}^{T}\text{diag}([\sum_{j=1}^{k}(M_{1}-M_{0})_{ij}])X_{C})\bigg{)}$
(157)
$\displaystyle\frac{\partial^{2}}{\partial\varepsilon}f_{1}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=\mu_{p}\bigg{(}\text{tr}(X_{C}^{T}(M_{1}-M_{0})X_{C})-\text{tr}(X_{C}^{T}\text{diag}([\sum_{j=1}^{k}(M_{1}-M_{0})_{ij}])X_{C})\bigg{)}$
(158)
$\displaystyle\frac{\partial^{2}}{\partial\epsilon^{2}}f_{1}(M_{0}+\varepsilon(M_{1}-M_{0})$
$\displaystyle=0$ (159)

Finding the directional derivative for $f_{2}$,

$\displaystyle f_{2}(M_{0}+\varepsilon(M_{1}-M_{0}))$ (160)
$\displaystyle=\frac{\mu_{p}}{b}\text{tr}\bigg{(}M_{0}+\varepsilon(M_{1}-M_{0})$
$\displaystyle-\text{diag}(\sum_{j=1}^{k}\big{(}M_{0}+\varepsilon(M_{1}-M_{0})\big{)}_{ij})\bigg{)}-\frac{1}{b}+k-k\log
b$ (161)
$\displaystyle\frac{\partial}{\partial\varepsilon}f_{2}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=\frac{\mu_{p}}{b}\text{tr}(M_{1}-M_{0})-\frac{\mu_{p}}{b}\text{tr}\bigg{(}\text{diag}(\sum_{j=1}^{k}(M_{1}-M_{0})_{ij})\bigg{)}$
(162)
$\displaystyle\frac{\partial^{2}}{\partial\varepsilon^{2}}f_{2}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=0$ (163)

Finding the directional derivative for $f_{3}$,

$\displaystyle f_{3}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=\frac{\text{tr}(M_{0}+\varepsilon(M_{1}-M_{0}))}{\tilde{P_{0}}}-\frac{\sum_{i=1}^{k}(\sum_{j=1}^{k}(M_{0}+\varepsilon(M_{1}-M_{0}))_{ij})^{2}}{\tilde{P_{0}}^{2}}$
(164)
$\displaystyle=\frac{\text{tr}(M_{0}+\varepsilon(M_{1}-M_{0}))}{\tilde{P_{0}}}-\frac{\sum_{i=1}^{k}(\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}(M_{0}+\varepsilon(M_{1}-M_{0}))_{iau})^{2}}{\tilde{P_{0}}^{2}}$
(165)
$\displaystyle\frac{\partial}{\partial\varepsilon}f_{3}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=$
$\displaystyle\frac{\text{tr}(M_{1}-M_{0})}{\tilde{P_{0}}}+\sum_{i=1}^{k}$
$\displaystyle
2(\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}(M_{1}-M_{0})_{iau})\times\frac{\sum_{i=1}^{k}(\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}(M_{0}+\varepsilon(M_{1}-M_{0}))_{iau})^{2}}{\tilde{P_{0}}^{2}}$
(166)
$\displaystyle\frac{\partial^{2}}{\partial\varepsilon^{2}}f_{3}(M_{0}+\varepsilon(M_{1}-M_{0}))$
$\displaystyle=\frac{\sum_{i=1}^{k}2(\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}(M_{1}-M_{0})_{iau})^{2}}{\tilde{P_{0}}^{2}}$
(167)

Adding up these three,

$\frac{\partial^{2}}{\partial\varepsilon^{2}}F(M_{0}+\varepsilon(M_{1}-M_{0}),\mathbf{t_{0}}+\varepsilon(\mathbf{t_{1}}-\mathbf{t_{0}}))=\frac{\sum_{i=1}^{k}2(\sum_{asu}\tilde{\pi}_{s}P_{as}x_{u}(M_{1}-M_{0})_{iau})^{2}}{\tilde{P_{0}}^{2}}$

which is continuous in $(M_{1},\mathbf{t_{1}})$ for all
$(M_{0},\mathbf{t_{0}})$ in a neighborhood of $(W,\pi)$.

###  D.3 Required Condition c): Upper bound of first derivative

With $G(S)=F(H(S),h(S))$, $\frac{\partial
G((1-\varepsilon)\mathbb{D}+\varepsilon
S)}{\partial\varepsilon}|_{\varepsilon=0^{+}}<-C<0\ \forall\ \pi,P$

$G(S)=f_{1}(H(S))+f_{2}(H(S))+f_{3}(H(S))$

Let
$\bar{S}=((S-\mathbb{D})t)P(\mathbb{D}t)^{T}+(\mathbb{D}t)P(S-\mathbb{D}t)^{T}$

For $f_{1}$

$\displaystyle f_{1}(H((1-\varepsilon)\mathbb{D}+\varepsilon S))$
$\displaystyle=\mu_{p}\text{tr}(X_{C}^{T}[H((1-\varepsilon)\mathbb{D}+\varepsilon
S)-\text{diag}(\sum_{j=1}^{k}H((1-\varepsilon)\mathbb{D}+\varepsilon
S)_{ij})]X_{C})$ (168) $\displaystyle H(S)$ $\displaystyle=(Sx)P(Sx)^{T}$
(169) $\displaystyle H((1-\varepsilon)\mathbb{D}+\varepsilon S)$
$\displaystyle=((1-\varepsilon)\mathbb{D}x+\varepsilon
St)P((1-\varepsilon)\mathbb{D}x+\varepsilon St)^{T}$ (170)
$\displaystyle=(\mathbb{D}x+\varepsilon(S-\mathbb{D})x)P(\mathbb{D}x+\varepsilon(S-\mathbb{D})x)^{T}$
(171)
$\displaystyle=(\mathbb{D}x)P(\mathbb{D}x)^{T}+\varepsilon(\mathbb{D}x)P((S-\mathbb{D})x)^{T}$
(172) $\displaystyle\quad\
+\varepsilon((S-\mathbb{D})x)P(\mathbb{D}x)^{T}+\varepsilon^{2}((S-\mathbb{D})x)P((S-\mathbb{D})x)^{T}$
(173) $\displaystyle\text{Finally,
}f_{1}((1-\varepsilon)\mathbb{D}+\varepsilon S)$
$\displaystyle=\text{tr}\big{(}X_{C}^{T}(\mathbb{D}x)P(\mathbb{D}x)^{T}X_{C}\big{)}+\varepsilon\
\text{tr}\big{(}X_{C}^{T}(\mathbb{D}x)P((S-\mathbb{D})x)^{T}X_{C}\big{)}$
(174) $\displaystyle\quad\ +\varepsilon\
\text{tr}\big{(}X_{C}^{T}((S-\mathbb{D})x)P(\mathbb{D}x)^{T}X_{C}\big{)}+\varepsilon^{2}\
\text{tr}\big{(}X_{C}^{T}((S-\mathbb{D})x)P((S-\mathbb{D})x)^{T}X_{C}\big{)}$
(175) $\displaystyle\quad\
+\text{tr}\big{(}X_{C}^{T}\text{diag}([H(\mathbb{D})_{i}])X_{C}\big{)}+\varepsilon^{2}\
\text{tr}\big{(}X_{C}^{T}\text{diag}([H(S-\mathbb{D})_{i}])X_{C}\big{)}$ (176)
$\displaystyle\quad\ +\varepsilon\
\text{tr}\big{(}X_{C}^{T}\text{diag}([\big{(}((S-\mathbb{D})x)P(\mathbb{D}x)^{T}+(\mathbb{D}x)P(S-\mathbb{D}x)^{T}\big{)}_{i}])X_{C}\big{)}$
(177) $\displaystyle\text{Now,}\frac{\partial
f_{1}}{\partial\varepsilon}\bigg{|}_{\varepsilon=0^{+}}$
$\displaystyle=\text{tr}\big{(}X_{C}^{T}(\bar{S}-\text{diag}([\bar{S}_{i}]))X_{C}\big{)}$
(178)

For $f_{2}$,

$\displaystyle f_{2}(H((1-\varepsilon)\mathbb{D}+\varepsilon S))$
$\displaystyle=\frac{\mu_{p}}{b}\text{tr}\bigg{(}H((1-\varepsilon)\mathbb{D}+\varepsilon
S)-\text{diag}([H((1-\varepsilon)\mathbb{D}+\varepsilon
S)_{i}])\bigg{)}-\frac{1}{b}+k-k\log b$ (179) $\displaystyle\frac{\partial
f_{2}}{\partial\varepsilon}\bigg{|}_{\varepsilon=0^{+}}$
$\displaystyle=\frac{\mu_{p}}{b}\text{tr}\bigg{(}\bar{S}-\text{diag}([\bar{S}_{i}])\bigg{)}$
(180)

For $f_{3}$,

$\displaystyle f_{3}(H((1-\varepsilon)\mathbb{D}+\varepsilon S))$
$\displaystyle=\frac{1}{\tilde{P}_{0}}\text{tr}\bigg{(}(\mathbb{D}x)P(\mathbb{D}x)^{T}+\varepsilon(\mathbb{D}x)P((S-\mathbb{D})x)^{T}$
(181) $\displaystyle\qquad\qquad\
+\varepsilon((S-\mathbb{D})x)P(\mathbb{D}x)^{T}+\varepsilon^{2}((S-\mathbb{D})x)P((S-\mathbb{D})x)^{T}\bigg{)}$
(182)
$\displaystyle-\frac{1}{\tilde{P}_{0}^{2}}\sum_{i=1}^{k}\sum_{j=1}^{k}\bigg{(}\bigg{[}(\mathbb{D}x)P(\mathbb{D}x)^{T}+\varepsilon(\mathbb{D}x)P((S-\mathbb{D})x)^{T}$
(183) $\displaystyle\qquad\qquad\
+\varepsilon((S-\mathbb{D})x)P(\mathbb{D}x)^{T}+\varepsilon^{2}((S-\mathbb{D})x)P((S-\mathbb{D})x)\bigg{]}_{ij}\bigg{)}^{2}$
(184) $\displaystyle\frac{\partial
f_{3}}{\partial\varepsilon}\bigg{|}_{\varepsilon=0^{+}}$
$\displaystyle=\frac{1}{\tilde{P}_{0}}\text{tr}\bigg{(}\bar{S}\bigg{)}-\frac{2}{\tilde{P}_{0}}\sum_{i=1}^{k}\bigg{(}\sum_{j=1}^{k}\bigg{[}(\mathbb{D}x)P(\mathbb{D}x)^{T}\bigg{]}_{ij}\times\sum_{j=1}^{k}[\bar{S}]_{ij}\bigg{)}^{2}$
(185)

From here, the proof is followed as in Appendix of Zhao et al. (2012) (in
Proof of Theorem 3.1) which also borrows from Bickel & Chen (2009).

##  Appendix E Dataset Summaries and Metrics

Name | p ($|\textbf{V}|$) | n ($|X_{i}|$) | e ($|E|$) | k ($y$)
---|---|---|---|---
Cora | 2708 | 1433 | 5278 | 7
CiteSeer | 3327 | 3703 | 4614 | 6
PubMed | 19717 | 500 | 44325 | 3
Coauthor CS | 18333 | 6805 | 163788 | 15
Coauthor Physics | 34493 | 8415 | 495924 | 5
Amazon Photo | 7650 | 745 | 238162 | 8
Amazon PC | 13752 | 767 | 491722 | 10
ogbn-arxiv | 169343 | 128 | 1166243 | 40
Brazil | 131 | 0 | 1074 | 4
Europe | 399 | 0 | 5995 | 4
USA | 1190 | 0 | 13599 | 4
Table 3: Datasets summary.

Refer to Table 3 for the dataset summary.

Metrics. A pair of nodes are said to be in agreement if they belong to the
same class and are assigned to the same cluster, or they belong to different
classes and have been assigned different clusters. For a particular
partitioning, ARI is the fraction of agreeable nodes in the graph. Accuracy is
obtained by performing a maximum weight bipartite matching between clusters
and labels. NMI measures the normalized similarity between the clusters and
the labels, and is robust to class imbalances. Mutual Information between two
labellings $X$ and $Y$ of the same data is defined as
$MI(X,Y)=\sum_{i=1}^{|X|}\sum_{j=1}^{|Y|}\frac{|X_{i}\cap
Y_{i}|}{N}\log\frac{N|X_{i}\cap Y_{i}|}{|X_{i}||Y_{i}|}$ and it is scaled
between 0 to 1.

##  Appendix F Training Details

All experiments were run on an NVIDIA A100 GPU and Intel Xeon 2680 CPUs. We
are usually running 4-16 experiments together to utilize resources (for
example, in 40GB GPU memory, we can run 8 experiments on PubMed
simultaneously). Again, the memory costs are more than dominated by the
dataset. All experiments used the same environment running CentOS 7, Python
3.9.12, PyTorch 2.0, PyTorch Geometric 2.2.0.

##  Appendix G VGAE

In a VGAE, the encoder learns mean ($\mathbf{\mu}$) and variance
($\mathbf{\sigma}$):
$\mu=\text{GCN}_{\mu}(X,A)\text{and}\log\sigma=\text{GCN}_{\sigma}(X,A)$ By
using the reparameterization trick, we get the distribution of the latent
space as:
$q(Z|X,A)=\prod_{i=1}^{N}q(\mathbf{z_{i}}|X,A)=\prod_{i=1}^{N}\mathcal{N}(\mathbf{z_{i}}|\mu_{i},\textrm{diag}(\sigma_{i}^{2}))$
A common choice for decoder is inner-product of the latent space with itself
which giving us the reconstructed $\hat{A}$.
$p(\hat{A}|Z)=\prod_{i=1}^{p}\prod_{j=1}^{p}p(\hat{A}_{ij}|z_{i},z_{j}),\ \
\text{with}\ \ p(\hat{A}_{ij}=1|z_{i},z_{j})=\textrm{sigmoid}(z_{i}^{T}z_{j})$

##  Appendix H Plots of evolution of latent space for other datasets

Refer to Figure 5. We can see the clusters forming in the latent space of the
VGAEs. In the case of Q-VGAE, since a GCN is used on this space, it can learn
non-linearities and the latent space shows different structures (like a
starfish in CiteSeer). Moreoever, these structures have their geometric
centres at the origin and grow out from there. In contrast, for Q-GMM-VGAE,
since a GMM is being learnt over the latent space, the samples are encouraged
to be normally distributed in their independent clusters, all of which have
different means and comparable standard deviations. So, we see multiple
"blobs", which more or less follow a normal distribution. This plot
effectively shows why a GMM-VGAE is more expressive than a VGAE.

Figure 5: Plots of evolution of latent space for Q-VGAE and Q-GMM-VGAE methods
for CiteSeer, PubMed, Brazil (Air Traffic) and Europe (Air Traffic) datasets.

##  Appendix I Attributed SBM theory and results

We validate the robustness and sensitivity of proposed methods to variance in
the node features and graph structure. We are also generating features using a
multivariate mixture generative model such that the node features of each
block are sampled from normal distributions where the centers of clusters are
vertices of a hypercube.

(a) Adjacency Matrix

(b) Feature Covariance
Matrix
$k_{f}=k$

(c) Feature Covariance
Matrix
$k_{f}>k$

(d) Feature Covariance
Matrix
$k_{f}<k$

Figure 6: Visualization of the generated adjacency and feature covariance
matrices for the ADC-SBM

SBM. The Stochastic Block Model (SBM)(Nowicki & Snijders, 2001) is a
generative model for graphs that incorporates probabilistic relationships
between nodes based on their community assignments. In the basic SBM, a
network with $p$ nodes is divided into $k$ communities or blocks denoted by
$C_{i}$, where $i=1,2,\cdots,k$. The SBM defines a symmetric block probability
matrix $B$ with size $(k\times k)$, where each entry $B_{ij}$ represents the
probability of an edge between a node in community $C_{i}$ and a node in
community $C_{j}$. Diagonal entries of this matrix represents the
probabilities of intra-cluster edges. This matrix $B$ captures the intra- and
inter-community connections and is assumed to be constant. $P(i\leftrightarrow
j|C_{i}=a,C_{j}=b)=B_{ab}$ denotes the probability of an edge existing between
nodes $i$ and $j$ when node $i$ belongs to community $a$ and node $j$ belongs
to community $b$. Using these probabilities, the SBM generates a network by
independently sampling the presence or absence of an edge for each pair of
nodes based on their community assignments and the block probability matrix
$B$.

Degree Corrected SBM. DC-SBM(Karrer & Newman, 2011) takes an extra set of
parameters $\theta_{i}$ controlling the expected degree of vertex $i$. Now,
the probability of an edge between two nodes (using the same notation as
above) becomes $\theta_{i}\theta_{j}B_{ab}$. This was introduced to handle the
heterogeneity of real-world graphs.

ADC-SBM Generation. We make use of the graph_tool library to generate the DC-
SBM adjacency matrix, with $p=1000,k=4$. To generate the $B$ matrix, we follow
the procedure in (Tsitsulin et al., 2023), by taking expected degree for each
node $d=20$ and expected sub-degree $d_{out}=2$. This gives us $B$ as:

$\begin{bmatrix}18&2&2&2\\\ 2&18&2&2\\\ 2&2&18&2\\\ 2&2&2&18\\\ \end{bmatrix}$

Also, $\theta$ is generated by sampling a power-law distribution with exponent
$\alpha=2$. We constrain the generated vector to $d_{min}=2$ and $d_{max}=4$.

To generate features, we use the make_classification function in the sklearn
library. We generate a $128$-dimensional feature vector for each node, with no
redundant channels. These belong to $k_{f}$ groups, where $k_{f}$ might not be
equal to $k$. We test three scenarios: a) matched clusters ($k_{f}=k$) b)
nested features ($k_{f}>k$) c) grouped features ($k_{f}<k$) as visualized in
Figure 6. Note that for better visualization, class_sep was increased to 5
(however, the results are given with a value of 1, which is a harder problem).

Additionally, we also consider both cases, with and without the coarsening
constraint term.

Results. Our objective is able to completely recover the ground truth labels
(NMI/ARI/ACC = 1) under all the specified conditions.

##  Appendix J Results on very large datasets

The results are presented in Table 4.

CoauthorCS CoauthorPhysics AmazonPhoto AmazonPC ogbn-arxiv Method ACC
$\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC $\uparrow$ NMI $\uparrow$ ARI
$\uparrow$ ACC $\uparrow$ NMI $\uparrow$ ARI $\uparrow$ ACC $\uparrow$ NMI
$\uparrow$ ARI $\uparrow$ ACC $\uparrow$ NMI $\uparrow$ ARI $\uparrow$ FGC
69.6 70.4 61.5 69.9 60.9 49.5 44.9 38.3 22.5 46.8 36.2 23.3 24.1 8.5 9.1
Q-MAGC (Ours) 70.2 76.4 60.2 75.3 67.2 66.1 70.4 66.6 58.6 62.4 51 31.1 35.8
24.4 15.6 Q-GCN (Ours) 85.4 79.6 79.7 85.2 72 81.6 66.3 57.6 48.3 56.7 42.4
28.8 34.4 27.1 19.7 Q-VGAE (Ours) 85.6 79.9 81.6 86.7 69 77.7 69.0 59.4 49.0
62.3 45.7 47.2 39.5 30.4 24.7 Q-GMM-VGAE (Ours) 70.1 72.5 61.6 83.1 71.5 76.9
76.8 67.1 58.3 55.5 56.4 40 OOM OOM OOM DMoN 68.8 69.1 57.5 45.4 56.7 50.3
61.0 63.3 55.4 45.4 49.3 47.0 25.0 35.6 12.7

Table 4: Comparison of methods on large attributed datasets.

##  Appendix K Implementation

The implementations for all the experiments can be found at
https://anonymous.4open.science/r/MAGC-8880/.

We have extensively used the PyTorch(Paszke et al., 2019) and PyTorch
Geometric(Fey & Lenssen, 2019) libraries in our implementations and would like
to thank the authors and developers.

##  Appendix L Explanation on why VAE manifolds are curved

Embedded manifolds obtained from VGAEs are curved and must be flattened before
any clustering algorithms using Euclidean distance are applied.

The latent space of a VAE is not constrained to be Euclidean. Connor et al.
(2021) point out that the variational posterior is selected to be a
multivariate Gaussian, and that the prior is modeled as a zero-mean isotropic
normal distribution which encourages grouping of latent points around the
origin. Works such as (Chen et al., 2020; Bogdanov & Shchur, 2021; Arvanitidis
et al., 2018) make the VAE latent space to be Euclidean/Hyperbolic/Riemannian,
and show good visuals. Moreover, it can be observed in our own work (Figure 3a
and Appendix Figure 5a) that the latent manifolds are curved and so, are not
suitable for conventional methods such as k-means clustering, which need
Euclidean distance.

##  Appendix M Complexities of some graph clustering methods

Some GCN _-based_ clustering methods:

  * •

AGC(Zhang et al., 2019) \- $\mathcal{O}(p^{2}nt+ent^{2})$
where $t$ is the number of iterations (within an epoch)

  * •

R-VGAE(Mrabah et al., 2022) \- $\mathcal{O}(pk^{2}n+(p(n+k)+e(p+k))$

  * •

S3GC(Devvrit et al., 2022) \- $\mathcal{O}(pn^{2}s)$
where $s$ is the average degree

  * •

HSAN(Liu et al., 2023b) \- $\mathcal{O}(pBn)$
they state it as $\mathcal{O}(B^{2}d)$ but that is only for 1 batch of size
$B$ and not the whole epoch

  * •

VGAECD-OPT(Choong et al., 2020) \- $\mathcal{O}(p^{2}nD^{L})$
where $D$ is the size of graph filter, $l$ is the number of linear layers

##  Appendix N Evolution of different loss terms throughout training

Figure 7 \- Each separate series has been normalized by its absolute minimum
value to see convergence behavior on the same graph easily. Every series is
decreasing/converging (except gamma, which represents sparsity regularization
and remains almost constant). Thus, we can be assured that no terms are
counteracting and hurting the performance. The legend is provided in the graph
itself. This plot is on the Cora dataset.

Figure 7: Evolution of the different loss terms throughout training, denoted
by their weight parameters. Also the term X_t^T theta_C X_t term is the
smoothness term tr$(X_{C}^{T}C^{T}\Theta CX_{C})$