# Extending bootstrap AMG for clustering of
attributed graphs

Pasqua D’Ambra  ,  Panayot S. Vassilevski  and  Luisa Cutillo
<EMAIL_ADDRESS>Institute for Applied Computing (IAC), CNR, Via P.
Castellino, 111, 80131, Napoli, Italy<EMAIL_ADDRESS>Fariborz Maseeh
Department of Mathematics and Statistics, Portland State University, Portland,
OR 97207, USA<EMAIL_ADDRESS>School of Mathematics, University of
Leeds, Leeds, LS2 9JT, UK

###### Abstract.

In this paper we propose a new approach to detect clusters in undirected
graphs with attributed vertices. We incorporate structural and attribute
similarities between the vertices in an augmented graph by creating additional
vertices and edges as proposed in [1, 2]. The augmented graph is then embedded
in a Euclidean space associated to its Laplacian and we cluster vertices via a
modified K-means algorithm, using a new vector-valued distance in the
embedding space. Main novelty of our method, which can be classified as an
early fusion method, i.e., a method in which additional information on
vertices are fused to the structure information before applying clustering, is
the interpretation of attributes as new realizations of graph vertices, which
can be dealt with as coordinate vectors in a related Euclidean space. This
allows us to extend a scalable generalized spectral clustering procedure which
substitutes graph Laplacian eigenvectors with some vectors, named
algebraically smooth vectors, obtained by a linear-time complexity Algebraic
MultiGrid (AMG) method.

We discuss the performance of our proposed clustering method by comparison
with recent literature approaches and public available results. Extensive
experiments on different types of synthetic datasets and real-world attributed
graphs show that our new algorithm, embedding attributes information in the
clustering, outperforms structure-only-based methods, when the attributed
network has an ambiguous structure. Furthermore, our new method largely
outperforms the method which originally proposed the graph augmentation,
showing that our embedding strategy and vector-valued distance are very
effective in taking advantages from the augmented-graph representation.

###### Key words and phrases:

attributed graphs, clustering, graph augmentation, bootstrap AMG

###### 1991 Mathematics Subject Classification:

05C50, 05C70, 65M55, 68T10

##  1\. Introduction

Graph clustering is a main task in network analysis [3, 4, 5]. In many real-
world networks, such as social and biomedical networks, in addition to the
graph structure some information related to possible vertex features are often
available; such additional information, often referred to as metadata, can be
used to improve the accuracy of clustering algorithms. Many clustering methods
have been extended to incorporate features information in the vertex
similarities detection so that the clusters are not only densely connected but
also very homogeneous in their attributes, as discussed in section 2. This
paper describes an overall approach that enables clustering of networks with
categorical vertex attributes. Our method builds on a recently proposed
generalized spectral method [6]. In the following we introduce the general
settings and briefly summarize our approach.

Given a set of vertices $V$ and a set of $m$ categorical labels
$L=\\{l_{1},\ldots,l_{m}\\}$, a vertex-attributed graph is the 4-tuple
$G=(V,E,{\mathcal{W}},L)$, where $E\subset V\times V$ is the edge set and
${\mathcal{W}}=(w_{ij})_{(i,j)\in E}$ is the (positive) edge weights matrix.
Each vertex $v_{i}\in V$ is associated with an attribute vector
$(l_{1}(v_{i}),\ldots,l_{m}(v_{i}))$, where $l_{j}(v_{i})$ is the attribute
value of vertex $v_{i}$ on attribute $l_{j}$. The cardinality or size of $G$
is the total number of vertices $|V|=n$.

In [6] we formulated the clustering of an undirected graph, as a partition
$V_{1},\ldots,V_{K}$, in $K$ nonempty sets, so that $V_{k}\subset V$,
$V_{k}\cap V_{h}=\emptyset,\ \forall k\neq h$ and $V_{1}\cup\ldots\cup
V_{K}=V$, which minimizes the following functional:

(1.1)
$RatioCut(V_{1},\ldots,V_{K})=\frac{1}{2}\sum_{k=1}^{K}\frac{W(V_{k},\overline{V_{k}})}{|V_{k}|},$

where $W(V_{k},\overline{V_{k}})=\sum_{i\in
V_{k},j\in\overline{V_{k}}}w_{ij}$, and $\overline{V_{k}}$ is the complement
of $V_{k}$ in $V$. We observe that a graph partition which minimizes (1.1)
corresponds to a partition which minimizes the weight of the edges between two
different sets and maximizes the number of the vertices within a set. The
solution of (a relaxed version of) this problem involves the eigenvectors of
the (weighted) graph Laplacian that are used to define embedding of the graph
into a Euclidean space of low dimension and use geometric means to construct
clusters [7]. In [6], we have generalized this concept by looking at the
spectrum of a generalized eigenvalue problem involving the (weighted) graph
Laplacian and an adaptive Algebraic MultiGrid (AMG) operator (linear solver)
which provided a natural way to select the dimension of the embedding
Euclidean space. Our main motivation was to design a linear complexity
algorithm which can be more scalable than standard spectral clustering
approaches which require computation of graph Laplacian eigenvectors. The
method demonstrated a very good ability to detect clusters in graphs with
medium and high modularities, i.e., when a block structure is well defined.

In this work, we extend the method to incorporate additional information on
graph vertices, so that the ability to detect clusters in graphs with a low
modularity can be improved while preserving the linear complexity. More
precisely, the aim is to re-formulate the method, to obtain partitions where
the vertices within one set have similar attribute values. In order to achieve
this purpose, we introduce an approach that starts with the augmentation of
the graph $G$ with new vertices representing the attribute values and new
edges connecting original vertices to new vertices representing their
attribute values. After this augmentation, we solve the minimization problem
associated to the relaxed form of the functional (1.1) (as in [6]) for the
augmented graph, this creates a low dimensional embedding space associated
with the Laplacian of the augmented graph. A modified K-means algorithm is
then applied in the embedding space, where a new distance measure between two
original vertices is defined.

The main original contribution and results of this paper can be summarised as
follows:

  * •

We propose a new generalized spectral clustering method to detect community in
attributed graphs; it belongs to the class of early fusion methods, where
attributes are fused to the structure information before applying the
embedding in a Euclidean space associated to some generalized eigenvectors of
the Laplacian graph. The new method allows us to extend the embedding
procedure introduced in [6], showing flexibility and robustness of the
procedure itself.

  * •

We introduce a new vector-valued distance in the embedding space, which has
its rationale into a block representation of the added (attribute) vertices
with their structure vertex origin, since they represent a particular feature
of their origin; then we extended the classical K-means by using the new
distance.

  * •

Extensive experiments on different types of synthetic datasets and real-world
attributed graph datasets widely used in the literature show that our new
algorithm employing attributes information improves the ability of the
structure-only-based methods, when the attributed network has an ambiguous
structure, with significant advantage with respect to state-of-the-art
structure-based methods. Comparisons with results available in the literature
also demonstrate that our new method largely outperforms the method which
originally proposed the graph augmentation employed in this paper, showing
that our embedding strategy and vector-valued distance are very effective in
taking advantages from the augmented-graph representation.

  * •

We finally made available all the functionalities used in our algorithm in a
freely-available software framework to be used in clustering of attributed
graphs.

The paper is organized as follows. In Section 2 we present related work on
various different approaches to fuse structure and attribute and classify our
method with respect to well recognized current taxonomy. In Section 3 we
summarize main features of a clustering algorithm already proposed by us to
partition un-directed graphs and whose extension to attributed graphs is the
objective of this work. In Section 4 we provide a description of a strategy
which augments the original graph with new vertices and new edges depending on
the attribute values associated to the original graph vertices.

In Section 5 we describe the architecture of the proposed software framework
implementing all the functionalities needed to apply the method and analyse
its asymptotic computational complexity. In Section 6 we present the
experimental setting and discuss results on synthetic graphs and on some real-
world data sets. Section 7 contains concluding remarks and possible future
works.

##  2\. Related Work

The growing need for embedding vertex attributes in graph clustering generates
a large amount of research methods aiming to fuse information on graph
structure and vertex features. Some interesting review are in [8, 9]. An
important issue is related to the correlation between vertex attributes and
graph structure. It is a general belief that the presence of structure-
attributes correlation improves the quality of the clustering, however some
recent works seem to refuse this assumption at least for some type of fusion
methods [10]. This discussion is beyond the scope of our work, and we focus on
analyzing the state-of-the-art techniques able to add vertex attributes in the
cluster detection. In [8] the authors classify the proposed fusion methods
into different categories depending on the adopted clustering approach, and
indeed many authors, driven by the clustering method, agree in classifying the
methods in two main classes named distance-based and model-based [11, 12, 13].
The aim of distance-based methods is to define a form of unified distance
which combine structural and attribute information, and then leverage on
existing clustering methods, such as K-means, functional optimization,
spectral clustering, to cluster attributed graphs based on the unified
distance. The unified distance can be defined in an implicit way, e.g., by
introducing weights on nodes [14], or explicitly, e.g., by distance functions
[15, 2, 1]. Model-based methods mainly rely on statistical inference and setup
a generative model which well fits both the structural and vertex attributes
information [11, 12, 16].

In [9] we can find an interesting prior classification depending on when
structure and attributes are fused with respect to the clustering algorithm.
Although the author focuses on some aspects of social networks, his review is
very general and exhaustive. Three overall classes of methods are considered:
early fusion methods, if structure information and attributes are fused before
applying a clustering algorithm, simultaneous fusion methods that fuse
structure and attributes simultaneously while applying the clustering
algorithm and finally, late fusion methods that apply clustering to the graph
structure and to the attribute labels separately and finally fuse the obtained
groups. The definition of these three classes is useful to understand the
level of influence of the attribute information on the method, and it
clarifies that methods belonging to the early fusion and late fusion classes
have often the advantage to reuse approaches that already demonstrated their
effectiveness in clustering of graphs with no use of additional information on
vertices. The methods are also classified in terms of the methodology applied
to find clusters. For the detailed classification we refer the reader to [9],
which also includes tables of main methods belonging to each class, and an
interesting ranking of the most influential methods in each class, based on
their PageRank values. On the base of the above taxonomy, our approach belongs
to the class of early fusion methods, and going down to the technique-related
subclasses, we can classify it in the subclass of node-augmented graph-based.
Main methods belonging to this class are the SA-cluster [2, 1] and its
improved versions. We will give details on the above methods in section 4. We
also observe that, our approach actually combines different methodologies
intersection of three subclasses, including also the sub-classes of embedding-
based and distance-based methods, as it is described in the following
sections.

##  3\. Embedding in the space of algebraically smooth vectors for graph
Laplacian

In [6] we proposed a new method for efficient embedding of undirected graphs
in a low-dimensional Euclidean space spanned by an accurate approximation of
the smallest part of the spectrum of a linear operator associated to the graph
Laplacian. The method can be viewed as an extension of the classical spectral
clustering which employs a standard eigenvalue problem for graph Laplacian.

Let ${\mathcal{D}}=diag(d_{i})_{i=1}^{n}$, with $d_{i}=\sum_{j=0}^{n}w_{ij}$,
be the diagonal matrix of weighted vertex degrees of $G$, the graph Laplacian
of $G$ is the matrix
${\mathcal{L}}={\mathcal{D}}-{\mathcal{W}}\in{\mathcal{R}}^{n\times n}$. We
observe that ${\mathcal{L}}$ is a symmetric, positive semi-definite M-matrix,
and its spectrum is a main tool for studying the graph [17]. In particular, in
spectral clustering, a graph is embedded in the Euclidean space spanned by the
first $K$ eigenvectors corresponding to the first $K$ smallest eigenvalues of
${\mathcal{L}}$; this embedding allows to apply standard spatial clustering
algorithms, e.g., the K-means optimization algorithm, giving good
approximations of optimal partitioning for well-clustered graphs [18].

As low-dimensional space for graph embedding, we use the space spanned by some
vectors, referred to as algebraically smooth vectors generated by an iteration
process. We use the bootstrap process coming from an Algebraic MultiGrid (AMG)
operator associated to the graph Laplacian. Note that we always consider the
graph Laplacian of a connected graph; in the case of graphs with more than one
connected component, we apply our method to each of its connected components.
Furthermore, we also eliminate singularity of the graph Laplacian by a rank-1
update which gives symmetric positive definite (s.p.d.) graph Laplacian matrix
${\mathcal{L}}_{S}={\mathcal{L}}+\lambda\;e\cdot e^{T}$, where $e$ is a vector
of dimension $n$ having unit (nonzero) components only for a single pair of
indices $i$ and $j$ corresponding to an arbitrarily chosen edge $(i,j)\in E$.
In the following, we briefly summarize the theoretical basis of our method and
its formulation.

In the theory of AMG [19], an $\epsilon-$smooth vector of the matrix
${\mathcal{L}}_{S}$ is defined as follows:

######  Definition 3.1.

Let ${\mathcal{M}}$ be a s.p.d. smoothing operator, i.e., the iteration matrix
of a relaxation method convergent in the ${\mathcal{L}}_{S}$-norm111The
${\mathcal{A}}-$norm of the vector $u$, where ${\mathcal{A}}$ is a s.p.d.
matrix is defined as $\|u\|_{{\mathcal{A}}}=\sqrt{u^{T}{\mathcal{A}}u}$., and
given $\epsilon\in(0,1)$, we say that the vector $u$ is algebraically
$\epsilon$-smooth with respect to ${\mathcal{L}}_{S}$ if

$\|(I-{\mathcal{M}}^{-1}{\mathcal{L}}_{S})u\|^{2}_{{\mathcal{L}}_{S}}\geq(1-\epsilon)\|u\|^{2}_{{\mathcal{L}}_{S}}.$

We note that
$\|(I-{\mathcal{M}}^{-1}{\mathcal{L}}_{S})u\|_{{\mathcal{L}}_{S}}\leq\|u\|_{{\mathcal{L}}_{S}}$
(since ${\mathcal{M}}$ is a convergent smoother). That is, $u$ is in the near-
nullspace of ${\mathcal{M}}^{-1}{\mathcal{L}}_{S}$ in the sense that
${\mathcal{M}}^{-1}{\mathcal{L}}_{S}u$ has a much smaller norm than $u$.

Such a vector can be obtained by performing a few iterations of the smoother
${\mathcal{M}}$, i.e.

(3.1) $u=u^{t}=(I-{\mathcal{M}}^{-1}{\mathcal{L}}_{S})u^{t-1},$

starting from a random vector $u^{0}$ for a sufficiently large $t$.

The above concept was generalized in [6] in the sense that we used an operator
${\mathcal{B}}$ (in fact a sequence of operators that gradually approximate
${\mathcal{L}}_{S}$) in place of ${\mathcal{M}}$ and considered the respective
algebraically smooth vectors.

Given the matrix ${\mathcal{L}}_{S}$, we first define a two-level method for
solving the linear system ${\mathcal{L}}_{S}x=b$, where
$x\in{\mathcal{R}}^{n}$ is the unknown vector and $b\in{\mathcal{R}}^{n}$ is a
given right-hand side term. Let $V_{c}\subset{\mathcal{R}}^{nc}$ be the space
spanned by a set of smooth vectors $\\{u_{r}\\}_{r}$ of ${\mathcal{L}}_{S}$
with respect to the operator ${\mathcal{M}}$. Let
$\\{\phi_{i}\\}_{i=1}^{n_{c}}$ provide a basis of $V_{c}$, and let
${\mathcal{P}}=[\phi_{1},\;\dots,\;\phi_{n_{c}}]$ be the interpolation matrix
mapping vectors from $V_{c}$ to $V\subset{\mathcal{R}}^{n}$ and
$({\mathcal{L}}_{S})_{c}={\mathcal{P}}^{T}{\mathcal{L}}_{S}{\mathcal{P}}$ the
corresponding coarse matrix. We define a two-level operator ${\mathcal{B}}$
corresponding to the iteration matrix

$I-{\mathcal{B}}^{-1}{\mathcal{L}}_{S}=(I-{\mathcal{M}}^{-T}{\mathcal{L}}_{S})(I-{\mathcal{P}}({\mathcal{L}}_{S})^{-1}_{c}{\mathcal{P}}^{T}{\mathcal{L}}_{S})(I-{\mathcal{M}}^{-1}{\mathcal{L}}_{S}),$

which leads to the formula

${\mathcal{B}}^{-1}=\overline{{\mathcal{M}}}^{-1}+(I-{\mathcal{M}}^{-T}{\mathcal{L}}_{S}){\mathcal{P}}({\mathcal{L}}_{S})^{-1}_{c}{\mathcal{P}}^{T}({\mathcal{I}}-{\mathcal{L}}_{S}{\mathcal{M}}^{-1}),$

where
$\overline{{\mathcal{M}}}={\mathcal{M}}({\mathcal{M}}+{\mathcal{M}}^{T}-{\mathcal{L}}_{S})^{-1}{\mathcal{M}}^{T}$
is the so-called symmetrized smoother, such that
${\mathcal{I}}-\overline{{\mathcal{M}}}^{-1}{\mathcal{L}}_{S}=({\mathcal{I}}-{\mathcal{M}}^{-T}{\mathcal{L}}_{S})({\mathcal{I}}-{\mathcal{M}}^{-1}{\mathcal{L}}_{S})$.

From the general theory of AMG methods [19] we can infer the following
spectral equivalence result for the operator ${\mathcal{B}}$:

$v^{T}{\mathcal{L}}_{S}v\leq v^{T}{\mathcal{B}}v\leq
C\;v^{T}{\mathcal{L}}_{S}v,\ \ \forall v\in V$

where the optimal constant $C$ is given by the formula:

(3.2) $C=\max\limits_{v\in V}\;\frac{\|v-\pi
v\|^{2}_{\mathcal{M}}}{\|v\|^{2}_{{\mathcal{L}}_{S}}}.$

Here $\pi:V\;\mapsto V_{c}\subset V$ is a projection in the
${\mathcal{M}}-$inner product, i.e.,
$\pi={\mathcal{P}}({\mathcal{P}}^{T}{\mathcal{M}}{\mathcal{P}})^{-1}{\mathcal{P}}^{T}{\mathcal{M}}$.
Furthermore, the following theorem is proved in [6]:

######  Theorem 3.1.

Consider the operator ${\mathcal{B}}$ defined by the two-level algorithm with
smoother ${\mathcal{M}}$, symmetrized smoother $\overline{{\mathcal{M}}}$, and
coarse space $V_{c}$. Let $\rho$ be the convergence rate of the iterative
process having ${\mathcal{B}}$ as iteration matrix. Consider the generalized
eigenvalue problem:

(3.3) ${\mathcal{L}}_{S}q_{k}=\lambda_{k}{\mathcal{M}}q_{k},$

where
$0\leq\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{s}<\lambda_{s+1}\leq\dots\leq\lambda_{\max}$,
and assume that the coarse space $V_{c}$ contains the first $s\geq 1$
eigenvectors $q_{k}$. Then, the following estimate holds for the spectral
equivalence constant $C$:

$C\leq\frac{1}{1-\rho}\simeq\frac{1}{\lambda_{s+1}}.$

The above result assures that if ${\mathcal{B}}$ is a method with a
sufficiently small convergence rate, then eigenvectors in (3.3) can
effectively be replaced (as argued in [6]) by the algebraically smooth vectors
generated by the iteration process based on ${\mathcal{B}}$. More
specifically, in [6] we proposed an iterative procedure which justified the
use, in place of the eigenvectors $q_{k}$, of a sequence of smooth vectors
with respect to ${\mathcal{L}}_{S}$. These vectors are obtained during a
bootstrap procedure whose final aim was to construct a composite AMG operator
${\mathcal{B}}$ leading to a sufficiently small spectrally equivalence
constant $C$. For sake of completeness, we summarize here main features of the
bootstrap procedure and the resulting graph embedding algorithm.

The bootstrap procedure builds a composite linear solver
$\overline{{\mathcal{B}}}$ obtained from a multiplicative composition of a
number of AMG operators ${\mathcal{B}}_{r}$. The operator
$\overline{{\mathcal{B}}}$ is defined from the error propagation matrix:

(3.4)
$I-\overline{{\mathcal{B}}}^{-1}{\mathcal{L}}_{S}=(I-{\mathcal{B}}^{-1}_{0}{\mathcal{L}}_{S})(I-{\mathcal{B}}^{-1}_{1}{\mathcal{L}}_{S})\dots\;(I-{\mathcal{B}}^{-1}_{r}{\mathcal{L}}_{S}),$

where each ${\mathcal{B}}_{r}$ is an AMG operator, e.g., a standard V-cycle
involving a simple smoother such as the weighted Jacobi or the Gauss-Seidel
method, and a hierarchy of coarse matrices exploiting coarsening based on
compatible weighted matching in the form studied in [20, 21, 22]. To construct
the operator ${\mathcal{B}}_{r}$, we first generate an algebraically smooth
vector $u_{r}$. The first smooth vector is constructed based on the smoother
${\mathcal{M}}$ for ${\mathcal{L}}_{S}$ as described in (3.1). At step $r>0$,
we apply few iterations based on the currently available composite solver
${\overline{\mathcal{B}}}$ consisting of the already built $r-1$ components;
i.e., $u=u^{t}=(I-\overline{{\mathcal{B}}}^{-1}{\mathcal{L}}_{S})u^{t-1}$,
starting from a random vector $u^{0}$. Once the vector $u_{r}$ is available,
the operator ${\mathcal{B}}_{r}$ is defined as an AMG cycle where the
construction of its hierarchy of coarse spaces is driven by the vector $u_{r}$
using coarsening based on compatible weighted matching [20, 21, 22].

The process is iterated several times until $M$ vectors
$\\{u_{r}\\}_{r=1}^{M}$ are computed; the number $M$ of computed vectors
depends on the desired convergence factor $\varrho$ that we impose on the
composite solver ${\mathcal{B}}_{BAMG}=\overline{{\mathcal{B}}}$ in (3.4) to
exhibit after $M$ composite steps. As a result, for a given $\varrho\in(0,1)$,
the number $M$ of components (equivalently the number of smooth vectors
associated with the components ${\mathcal{B}}_{r}$) are such that by
construction, we have:

$\|{\mathcal{I}}-{\mathcal{B}}^{-1}_{BAMG}{\mathcal{L}}_{S}\|_{{\mathcal{L}}_{S}}\leq\varrho.$

In order to generate a set of orthogonal vectors spanning the graph embedding
space, we finally apply a Singular Value Decomposition (SVD) to the computed
set of vectors $\\{u_{r}\\}^{M}_{r=1}$ and use the corresponding left-singular
vectors $\\{\overline{u}_{r}\\}_{r=1}^{n_{c}}$ as basis for the embedding
space. In addition to this set we add the vector of all ones, which is the
null vector of ${\mathcal{L}}$, and we also use this vector as the starting
one in the bootstrap process. After embedding, we apply any standard spatial
clustering procedure to identify possible vertex clusters of the original
graph. Our method of choice is the widely used K-means algorithm, as detailed
in [6] (see Algorithm 2).

Note that our procedure exploits increasingly accurate coarse spaces which
provide improved approximations to the generalized eigenvalue problem (3.3).
This is the case since the operator ${\mathcal{B}}$ changes at each bootstrap
step until we reach a final operator which has a sufficiently small spectral
equivalence constant (3.2). The achieved convergence tolerance determines how
many bootstrap steps to perform, and as a result gives a criterion how to
choose the dimension of the embedding space. For the more standard spectral
clustering it is not as clear how many eigenvectors to compute in advance in
order to get satisfactory clustering results. Therefore, our method exhibits a
desired property in practice and we view it as an advantage over the more
standard spectral clustering approaches.

In the following, we describe the extension of the method based on the
bootstrap AMG to undirected graphs with vertex attributes.

##  4\. Attribute-augmented Graphs

Let $G=(V,E,{\mathcal{W}},L)$ be a vertex-attributed graph, where $L$ is the
set of attribute values, i.e., the set of vectors
$(l_{1}(v_{i}),\ldots,l_{m}(v_{i}))$ with $m$ values associated to each vertex
$v_{i}\in V$, and let $Dom(l_{j})=\\{l_{j,1},\ldots,l_{j,n_{j}}\\}$ be the
$n_{j}$ distinct values of the attribute $l_{j}$. Algorithm 1 builds an
augmented graph with additional vertices $V_{a}$, edges $E_{a}$ and edge
weights $W_{a}$. It implements the method from [2, 1], and can be viewed also
as a coarsened version of an attribute-based graph disaggregation algorithm
[23].

Data: vertex-attributed graph $G=(V,E,{\mathcal{W}},L)$, $n=|V|$, $ne=|E|$,
$m$: number of label classes in $L$, $weight(\cdot)$: real function depending
on vertices and attribute values

Result: augmented graph $G_{A}^{H}=(V\cup V_{a},E\cup
E_{a},{\mathcal{W}}_{a})$, $n_{new}=|V\cup V_{a}|$, $ne_{new}=|E\cup E_{a}|$

for  _$j=1,\ldots,m$_ do

identify $Dom(l_{j})=\\{l_{j,1},\ldots,l_{j,n_{j}}\\}$ as the set of distinct
values of the label $l_{j}$;

end for

$E_{a}=\emptyset$, $V_{a}=\emptyset$, ${\mathcal{W}}_{a}={\mathcal{W}}$;

for  _$j=1,\ldots,m$_ do

for  _$k=1,\ldots,n_{j}$_ do

$V_{a}=V_{a}\cup\\{v_{j,k}\\}$;

for  _$i=1,\ldots,n$_ do

if  _$l_{j}(v_{i})==l_{j,k}$_ then

$E_{a}=E_{a}\cup\\{e_{new}=(v_{i},v_{j,k})\\}$;

${\mathcal{W}}_{a}(e_{new})=weight(e_{new})$;

end if

end for

end for

end for

$n_{new}=n+\sum_{j=1}^{m}n_{j}$;

$ne_{new}=ne+nm$;

ALGORITHM 1 Graph augmentation based on vertex attributes.

In summary, Algorithm 1 builds the augmented graph defined as follows:

######  Definition 4.1.

Let $G=(V,E,{\mathcal{W}},L)$ be a vertex-attributed graph, where $L$ is the
set of attribute values, i.e., the set of vectors
$(l_{1}(v_{i}),\ldots,l_{m}(v_{i}))$ with $m$ values associated to each vertex
$v_{i}\in V$, and let $Dom(l_{j})=\\{l_{j,1},\ldots,l_{j,n_{j}}\\}$ be the
$n_{j}$ distinct values of the attribute $l_{j}$. The augmented graph
$G_{A}^{H}=(V\cup V_{a},E\cup E_{a},{\mathcal{W}}_{a})$, constructed in
Algorithm 1 has its components defined as:

  1. (1)

$V_{a}=(v_{j,k})$ is the set of vertices representing the distinct attribute
values $(l_{j,k})_{j=1,k=1}^{m,n_{j}}$ of the vertices in $G$;

  2. (2)

$E_{a}$ is the set including the edges $(v_{i},v_{j,k})\ \forall j=1,\ldots,m$
iff $l_{j}(v_{i})=v_{j,k}$, which connect the vertex $v_{i}$ with his
attribute value $v_{j,k}$;

  3. (3)

${\mathcal{W}}_{a}$ is a real matrix including possible weights for the edges
in $E\cup E_{a}$.

We note that the cardinality of the augmented graph $G^{H}_{A}$ is equal to
$n_{new}=n+\sum_{j=1}^{m}n_{j}$, where $n_{j}\leq n$ is the number of
different values of the attribute $l_{j}$, and the number of edges is
$ne_{new}=ne+nm$, where $ne$ is the number of edges in $G$. We also observe
that in our implementation of the method and in the experiments discussed in
section 6 we consider unweighted graphs, where matrix ${\mathcal{W}}$ and
${\mathcal{W}}_{a}$ corresponds to the adjacency matrix of the original graph
and of the augmented graph, respectively. On the other hand, the extension of
the code to the weighted graphs is straightforward, while the impact of
possible weights on the edges could have large impact on the final clusterings
for any employed method.

Figure 1. Original (left) and Augmented Graph (right).

In Fig. 1 we show a small undirected graph with $3$ different labels and the
corresponding augmented graph. The original graph has the vertex set
$V=\\{v_{1}=1,v_{2}=2,v_{3}=3,v_{4}=4\\}$ and the edge set
$E=\\{(1,2),(1,3),(1,4),(2,4),(2,3)\\}$, then $|V|=n=4$ and $|E|=ne=5$; the
vertices are characterized by the labels defined in Table 1.

vertex | labels
---|---
| $l_{1}$ | $l_{2}$ | $l_{3}$
1 | M | R | C
2 | F | D | P
3 | F | I | J
4 | M | D | C
Table 1. Sample Graph: Vertex Attributes

In this case, we have $m=3$, $Dom(l_{1})=\\{M,F\\}$, $Dom(l_{2})=\\{R,D,I\\}$,
$Dom(l_{3})=\\{C,P,J\\}$, $n_{1}=2$ and $n_{2}=n_{3}=3$. Algorithm 1 builds
the set of $8$ new vertices, enumerated with consecutive integers from $5$ to
$12$. We have
$V_{a}=\\{v_{1,1}=M=5,v_{1,2}=F=6,v_{2,1}=R=7,v_{2,2}=D=8,v_{2,3}=I=9,v_{3,1}=C=10,v_{3,2}=P=11,v_{3,3}=J=12\\}$.
The set of new edges reads:
$E_{a}=\\{(1,5),(1,7),(1,10),(2,6),(2,8),(2,11),(3,6),(3,9),(3,12),(4,5),(4,8),(4,10)\\}$.
Finally, the augmented graph $G^{H}_{A}$ has a total of $n_{new}=12$ vertices
and a total of $ne_{new}=17$ edges. Both, the original graph and the augmented
graphs have no edge weights.

In the augmented graph (see Definition 4.1) two original vertices $v_{i_{1}}$
and $v_{i_{2}}$, referred to as structure vertices, can be close either if
they are connected through other structure vertices, or if they share
attribute vertices as neighbors, or both. In [2, 1] the authors have defined a
unified distance measure on the augmented graph using random walks, which
estimates the pairwise vertex closeness in the graph through both structure
and attribute edges. Using their distance measure, they apply a modified
K-means algorithm constituting the so-called Structural/Attribute clustering
algorithm (SA-cluster).

In this paper, we follow a different approach which can be seen as a
generalization of a spectral clustering algorithm extending the method
proposed in [6] to deal with vertex-attributed graphs. In details, after
augmentations we compute the s.p.d modified graph Laplacian
${\mathcal{L}}_{s}$ of the augmented graphs $G^{H}_{A}$ and embed the graphs
in the subspace spanned by the smooth vectors of ${\mathcal{L}}_{S}$
associated to a composite AMG operator $B_{BAMG}$ as in (3). These smooth
vectors, having $n_{new}$ entries are used to represent the original $n$
structure vertices in terms of a set of $1+m$ scalar components, where the
first component is related to the original graph structure information and the
resulting $m$ are related to the vertex attributes. In this way, each original
(structure) vertex is associated to $m$ additional (attribute) vertices, and
can be viewed as a node of $m+1$ disaggregated vertices (in analogy to [23]).
Another analogy can be seen with the representation of a vector field in a
physical (2D or 3D) space domain where the value of the field in a space point
is represented by a node vector of (2 or 3) scalar components. The latter
point of view was our motivation to define in section 4.1 the distance metric
(4.1), which is another key ingredient of our method.

Before focusing on our approach to identify partitions in the structure
vertices after embedding of augmented graph, we point out that there is also
another potential feature in exploiting the augmented graph with the goal to
analyze important properties of the original vertex-attributed graph, that is
the detection of overlapped clusters. To illustrate the idea, consider our
example in Fig. 1. We can use any existing clustering technique applied to the
augmented graph, e.g., the method we studied in [6] or any other available
method. Assume that one of these methods partitioned the augmented graph into
the following two clusters:

${\mathcal{A}}^{a}_{1}=\\{1,4,5,7,8,10\\}\text{ and
}{\mathcal{A}}^{a}_{2}=\\{2,3,6,9,11,12\\}.$

If we restrict these sets to the original (structure) vertices, we get the
following two clusters:

${\mathcal{A}}_{1}=\\{1,4\\}\text{ and }{\mathcal{A}}_{2}=\\{2,3\\},$

which provides one possible way to construct clusters for vertex-attributed
graphs. In the present paper, we cluster the original graph by embedding the
augmented one into an Euclidean space of some dimension and then define
distance between two structure (original) vertices only, so that a new
distance-based K-means is applied to them. We observe that this approach is
computationally less expensive than applying classical K-means to all vertices
of the augmented graph. However, assume that we did cluster the augmented
graph leading to the clusters ${\mathcal{A}}^{a}_{1}$ and
${\mathcal{A}}^{a}_{2}$. It allows us to additionally exploit the attributed
vertices as explained in the following. If we replace the new (attribute)
vertices with their related structure vertices in the sets
${\mathcal{A}}^{a}_{1}$ and ${\mathcal{A}}^{a}_{2}$, i.e., $5\mapsto 1,4$;
$6\mapsto 2,3$; $7\mapsto 1$; $8\mapsto 2,4$; $9\mapsto 3$; $10\mapsto 1,4$;
$11\mapsto 2$, and $12\mapsto 3$, we end up with the following new
partitioning of the original structure vertices:

${\widetilde{\mathcal{A}}}_{1}=\\{1,2,4\\}\text{ and
}{\widetilde{\mathcal{A}}}_{2}=\\{2,3,4\\}.$

We notice that we identified two overlapping clusters. Therefore, the
augmentation process, if fully utilized, offers a potential to create
overlapping subgroups of structure (original) vertices.

###  4.1. Vertex block-coordinates and new vector-valued Euclidean distance
measure

####  4.1.1. Vertex block-coordinates

Starting from the augmented graph $G_{A}^{H}$ associated with the attributed
graph $G$, let
$\overline{u}_{r}=(\overline{u}^{r}_{1},\ldots,\overline{u}^{r}_{n},\ldots,\overline{u}^{r}_{n_{new}})\in{\mathcal{R}}^{n_{new}}$
be a smooth vector computed for the augmented graph Laplacian of $G_{A}^{H}$.
We recall that the vertices of the augmented graph come with double
subscripts; namely, structure vertex $v_{i}$ is indexed as $v_{i,1}$ whereas
all its attribute neighbors $v_{j,k}$ are indexed such that for attribute
$j=1,2,...,m$ we have $l_{j}(v_{i})=v_{j,k}$ (the value from the domain of
$l_{j}$). According to this convention (applying it to the indices of the
smooth vectors), we can partition the values of any smooth vector in blocks,
where each block will correspond to a structure vertex and its $m$ attribute
neighbors. That is, we can form block coordinate vectors
$v^{r}_{i}=(\overline{u}^{r}_{i,j})_{j=1}^{m+1}$ for each structure vertex
$v_{i}$ and each smooth vector $\overline{u}_{r}$ where, as already indicated,
we let $\overline{u}^{r}_{i,1}=\overline{u}^{r}_{i}$, and
$\overline{u}^{r}_{i,j+1}=\overline{u}^{r}_{j,k}$ if $l_{j}(v_{i})=v_{j,k}$,
$\forall j=1,\ldots,m$. In this way, the original graph vertices will have the
same coordinate $\overline{u}^{r}_{j,k}$ if they share the same attribute
value $v_{j,k}$ on a given attribute $l_{j}$, so that vertices sharing an
attribute value belong to the same hyperplane corresponding to that value.

In Fig. 2 we draw in a 3D space a picture of the coordinates representing the
$4$ original graph vertices obtained by using only $1$ smooth vector generated
by the bootstrap algorithm applied to the augmented graph Laplacian. We fix
the coordinate related to the structure information on $x$-axis and use a
combination of the other $3$ coordinates associated to the attributes as $y-$
and $z-$coordinates. We also draw the planes corresponding to the same
attribute values, so that vertices sharing attribute coordinates belong to the
same plane.

Figure 2. Coordinates representation for the sample graph: the $4$ coordinates
of each original graph vertex, corresponding to $1$ smooth vector, are
represented in 3D, fixing on $x-$axis the coordinate related to the structure
information and varying on $y-$axis and $z-axis$ the coordinates related to
different labels. Planes corresponding to the same attribute values are also
drawn.

####  4.1.2. Distance metric

After assigning to each structure vertex $v_{i}$ a block “augmented” smooth
vector component $v^{r}_{i}\in{{\mathcal{R}}}^{1+m}$, it is straightforward to
define a distance between any two structure vertices $v_{i_{1}}$ and
$v_{i_{2}}$ based on these $(1+m)-$vector coordinates.

Let $n_{c}$ be the number of smooth vectors obtained by the bootstrap AMG
applied to the augmented graph Laplacian ${\mathcal{L}}_{S}$ and let
$v_{i_{1}}=(\overline{u}^{r}_{i_{1},j})^{m+1,n_{c}}_{j=1,r=1}$ and
$v_{i_{2}}=(\overline{u}^{r}_{i_{2},j})^{m+1,n_{c}}_{j=1,r=1}$ be the block-
coordinates of the vertices, then we let

(4.1)
$dist(v_{i_{1}},v_{i_{2}})=\sqrt{\sum\limits^{n_{c}}_{r=1}\|\overline{u}^{r}_{i_{1}}-\overline{u}^{r}_{i_{2}}\|_{2}^{2}},$

where $\|\overline{u}^{r}_{i_{1}}-\overline{u}^{r}_{i_{2}}\|_{2}^{2}$ is the
square of the usual Euclidean distance between the block components
$(\overline{u}^{r}_{i_{1},j})^{m+1}_{j=1}$ and
$(\overline{u}^{r}_{i_{2},j})^{m+1}_{j=1}$.

The function in (4.1) is clearly a metric (due to the standard Cauchy-Schwarz
inequality for vectors).

Note that if we did not have labels, i.e., if $m=0$, the above expression
reduces to the standard definition of distance once we have embedded the graph
in a Euclidean space. For $m>0$, our motivation is to block the added vertices
with their structure vertex origin since they represent a particular feature
of their origin.

Once we have a distance measure defined, we are able to apply a K-means
clustering algorithm for the set of structure vertices.

##  5\. BCMatch4Graphs Software Framework

We implemented all the functionalities needed to apply our clustering
algorithms, both for simple graphs and for vertex attributed ones, extending
the software framework described in [21]. The extended code to deal with
graphs has been made available with the new name `BCMatch4Graphs` at
https://github.com/bootcmatch/BCMatch4Graphs.

Figure 3. Flowchart of a typical computation workflow based on BCMatch4Graphs.

In Fig. 3 we show a flowchart of a typical computation workflow supported by
the software framework whose architecture has a suitable layered structure. We
extended the basic layer of the modular BootCmatch framework with some
functionalities for manipulation of un-directed graphs, computation of related
matrix representations, including Laplacian matrix and rank-1 update to
eliminate constant vector from its kernel. In this layer we also included the
code for graph augmentation to be applied if vertex attributes are available,
as described in Algorithm 1. We note that a support function to extract from
each class of attributes only their different values, i.e., the set
$Dom(l_{j})$, per each class label $l_{j}$, is included in the augmentation
algorithm. All the functionalities for graph embedding are exposed from a
middle software layer implementing the setup and the application of the
operator in (3), to which correspond a number of algebraically smooth vectors
spanning the graph embedding space. After smooth vectors computation, a
function to compute the SVD is applied to define a vector basis. Finally,
functions needed to define vertex coordinates and then compute the usual
Euclidean distance for simple graphs and the new vector-valued distance as
defined in (4.1) for vertex attributed graphs, are also included in the middle
layer. At the application layer, functionalities to apply the K-means
iterative procedure and evaluate clustering quality, i.e., compute clustering
modularity and conditional entropy with respect to the set of attributes, as
defined in (6.2) and (6.2), respectively, are implemented. Some user’s level
drivers to cluster un-directed connected graph, with and without vertex
attributes, are also available in the framework.

###  5.1. Computational complexity

We observe that the way we build and apply each component of the bootstrap AMG
operator in (3) defines the computational complexity of the graph embedding
procedure and then its convergence properties, i.e., the number of components
needed to reach a desired convergence rate. The embedding software layer in
`BCMatch4Graphs` supports different AMG cycles among the most widely used,
featured by different computational complexities also depending on the
algorithm applied to generate each single AMG operator. Furthermore,
symmetrized versions of the multiplicative composition of the built AMG
components can be applied, producing a more accurate composite AMG. The choice
of accuracy of each AMG operator and the way we compose the different built
operators influences the number of bootstrap steps, and then the number of
computed smooth vectors needed to reach the desired convergence rate
decreases. This automatic way to compute the dimension of the embedding vector
space is one of the most important advantages of our approach, and the large
flexibility in selecting the algorithmic parameters to setup and apply the
composite operator $\overline{{\mathcal{B}}}$ defines the final computational
complexity of the overall clustering procedure, which can be tuned depending
on the user’s needs. Some best practice choices, as detailed in section 6, are
available so that the final embedding procedure has a linear computational
complexity $\mathcal{O}(M\cdot ne_{new})$, where $ne_{new}$ is the number of
edges in the augmented graph and $M$ is the number of bootstrap steps (and
then of the computed smooth vectors) needed to reach the desired convergence
rate. Therefore, the final computational complexity of the clustering method
is $\mathcal{O}(M\cdot ne_{new}+M\cdot n_{new}\cdot nit)$, where
$\mathcal{O}(M\cdot n_{new}\cdot nit)$ is the complexity of the K-means
procedure if $nit$ cluster iterations are performed.

##  6\. Experimental Results

In the present section, we study the clustering quality of our method
described in section 4 on two types of synthetic graphs widely used as
benchmarks and on some real-life attributed graphs available from public
repositories. We also compare our algorithm with some existing clustering
methods, whose implementation is freely available, some of which also use
information on vertices attributes. More specifically, we consider the
following attribute-free methods:

  * •

Louvain method ($Lou$) is a greedy optimization method that attempts to
maximize the modularity of a partion of the graph by a multilevel algorithm
[24];

  * •

Fast Greedy method ($FG$) is also a greedy method which tries to maximize
modularity by a hierarchical agglomeration algorithm [25];

  * •

$BCMG$ is our original method, available in `BCMatch4Graphs` and based on the
bootstrap AMG for the Laplacian graph [6]. We point out that some default
algorithmic parameters to setup the method are applied, as discussed in [6]
and defined in the same sample of configuration files available at the
software repository. In order to reduce sensitivity of the K-means algorithm
w.r.t. the starting cluster centers, we run the k-means for $100$ times with
the same input data and choose the partition corresponding to the maximum
achieved value of the computed modularity. In the application of the bootstrap
AMG, we built $m$ smooth vectors, where $m=\min(40,\overline{K})$, with
$\overline{K}$ the number of AMG operators built to obtain a final
multiplicative operator of type (3.4) with a desired convergence rate less
than $\varrho=10^{-8}$. Each AMG operator is a symmetric V-cycle with one
forward/backward Gauss-Seidel sweep as pre/post-smoothing. A direct LU
factorization is used as coarsest solver.

Then we consider the following clustering methods that employ information on
vertices attributes:

  * •

$Newman$ is the method described in [16]. It uses a techniques of Bayesian
statistical inference to fit a suitable network model to the available data,
trying to find a possible partition correlated to the attributes;

  * •

$BAGC$ is a model-based approach to attributed graph clustering. It also
relies on a Bayesian probabilistic model which tries to fit both structural
and attribute aspects of networks [11];

  * •

$BCMAG$ is the new clustering method based on the augmented graph described in
this paper and implemented in `BCMatch4Graphs`. The bootstrap algorithm for
augmented graph embedding and the modified K-means are applied in the same
conditions of the above $BCMG$ method.

Note that $Lou$ and $FG$ are implemented in the $R$ package $igraph$. $BAGC$
is available on the GitHub repository [26]. Software implementing $Newman$ is
available as supplementary material of [16]. We observe that this software
does not allow the use of more than one class of attributes. For this reason
we could not show any results from $Newman$ method in Section 6.3 and in some
of the test cases in Section 6.4.

###  6.1. Evaluation Metrics

In order to analyze the obtained results, we employ standard metrics which
allow us to compare two different clusterings and then quantify the
performance of the clustering algorithms when the ground truth is given. Let
${\mathcal{C}}^{\prime}$ be the ground truth, and ${\mathcal{C}}$ be the
estimated partition obtained by a given clustering method, in the following we
define Normalized Mutual Information (introduced in [27]), entropy and gain.

Given the graph $G$ and two partitions ${\mathcal{C}}={C_{1},\ldots,C_{K}}$
and ${\mathcal{C}}^{\prime}={C^{\prime}_{1},\ldots,C^{\prime}_{K^{\prime}}}$,
with $K$ and $K^{\prime}$ non empty sets, respectively, we consider

  * •

$p(k)$ to be the probability of a vertex to belong to group $C_{k}$.

  * •

$p(k,k^{\prime})$ to be the probability that a vertex belongs to set $C_{k}$
from ${\mathcal{C}}$ and to the group $C^{\prime}_{k^{\prime}}$ from partition
${\mathcal{C}}^{\prime}$.

If we set $n_{k}=|{\mathcal{C}}_{k}|$ and
$n_{k,k^{\prime}}=|{\mathcal{C}}_{k}\cap{\mathcal{C}}_{k^{\prime}}|$, we can
define $p(k)=n_{k}/n$ and $p(k,k^{\prime})=n_{k,k^{\prime}}/n$.

$NMI$ is defined as:

$NMI({\mathcal{C}},{\mathcal{C}}^{\prime})=\frac{I({\mathcal{C}},{\mathcal{C}}^{\prime})}{H({\mathcal{C}})+H({\mathcal{C}}^{\prime})},$

where $H(C)$ is the entropy associated with partition ${\mathcal{C}}$:

(6.1) $H({\mathcal{C}})=-\sum_{k=1}^{K}p(k)\log(p(k)),$

and $I(C,C^{\prime})$ is the mutual information between ${\mathcal{C}}$ and
${\mathcal{C}}^{\prime}$, i.e., the information that one partition has about
the other:

$I({\mathcal{C}},{\mathcal{C}}^{\prime})=\sum_{k=1}^{K}\sum_{k^{\prime}=1}^{K^{\prime}}p(k,k^{\prime})\log\frac{p(k,k^{\prime})}{p(k)p(k^{\prime})}.$

Note that $NMI\in[0,1]$, and two partitions can be considered in good
agreement when $NMI\approx 1$.

Given the entropy $H({\mathcal{C}})$ and $H({\mathcal{C}}^{\prime})$ of two
clusterings, the cluster-specific entropy of ${\mathcal{C}}^{\prime}$ , that
is the conditional entropy of ${\mathcal{C}}^{\prime}$ with respect to cluster
$C_{k}$, is defined as:

$H({\mathcal{C}}^{\prime}|C_{k})=-\sum_{k^{\prime}=1}^{K^{\prime}}\frac{n_{k,k^{\prime}}}{n_{k}}\log\left(\frac{n_{k,k^{\prime}}}{n_{k}}\right),$

The conditional entropy of ${\mathcal{C}}^{\prime}$ given clustering
${\mathcal{C}}$ is then defined as the weighted sum:

(6.2)
$H({\mathcal{C}}^{\prime}|{\mathcal{C}})=\sum_{k=1}^{K}\frac{n_{k}}{n}H({\mathcal{C}}^{\prime}|C_{k})=-\sum_{k=1}^{K}\frac{n_{k}}{n}\sum_{k^{\prime}=1}^{K^{\prime}}\frac{n_{k,k^{\prime}}}{n_{k}}\log\left(\frac{n_{k,k^{\prime}}}{n_{k}}\right).$

Expression (6.2) suggests that the more cluster’s members are split into
different partitions, the higher the conditional entropy. If
${\mathcal{C}}^{\prime}$ represents the ground truth, for a perfect clustering
the conditional entropy value, which for the sake of brevity we refer to as
entropy in the rest of the paper, should be zero.

Suppose we want to measure if a certain partition reduces the overall entropy
and hence is more informative. The measure that we suggest to look at is
called information gain, or more simply, gain. Information gain is the
expected reduction in entropy caused by partitioning the vertices according to
a given partition ${\mathcal{C}}$. It is defined as the difference between the
absolute entropy of the clustering ${\mathcal{C}}^{\prime}$ and the
conditional entropy $H({\mathcal{C}}^{\prime}|{\mathcal{C}})$:

$IG({\mathcal{C}}^{\prime},{\mathcal{C}})=H({\mathcal{C}}^{\prime})-H({\mathcal{C}}^{\prime}|{\mathcal{C}})$

###  6.2. SBM Graphs

In order to assess the overall proposed methodology we started from a
Stochastic Block Model (SBM) [28] to generate random graphs. This type of
synthetic graphs are often used to identify specific model parameters that are
critical for the ability of an algorithm to find reliable partitions.

Assume there are $q$ sets of vertices in a graph, a SBM is specified by the
expected fraction of vertices in each set $n_{a},~{}1\leq a\leq q$, and by the
probability $p_{ab}$ of an edge between set $a$ and set $b$ $\forall
a,b\in[1,\ldots,q]$. These probabilities form a matrix $P$ that is commonly
referred to as affinity matrix. We focus on a simple well-known (cf., [29])
SBM model called planted graph. The planted model is a special SBM case in
which the values of the probability in the affinity matrix $P$ are equal to a
constant $p_{in}$ on the diagonal and equal to another constant $p_{out}$ off
the diagonal; thus, two vertices within the same set share an edge with
probability $p_{in}$, while two vertices in different sets share an edge with
probability $p_{out}$. In particular, if $p_{in}>p_{out}$ the model is called
assortative, while if $p_{in}<p_{out}$ the model is called dis-assortative. In
our simulations we assume an assortative structure, so that $p_{in}>p_{out}$.
We generate a SBM random graph $G$ of $n$ vertices partitioned into $q$ sets,
$C_{1},\ldots,C_{q}$, with adjacency matrix ${\mathcal{A}}_{ij}=1$ if there is
an edge between vertices $i$ and $j$, and $0$ otherwise; more specifically, we
consider Bernoulli trials for each potential edge with the probabilities given
by the affinity matrix, so that the corresponding adjacency matrix has entries
defined as follows:

${\mathcal{A}}_{rs}|~{}r\in C_{a}\>s\in C_{b}\sim Bern(p_{ab}),~{}\forall
r,s\in[1,\ldots,n]$

where for $a,b\in[1,\ldots,q]$:

$p_{ab}=\begin{cases}p_{in}&\text{if $a=b$ }\\\ p_{out}&\text{otherwise}\\\
\end{cases}$

Note that the probabilities are defined in such a way that the average vertex
degree is a given $c$. In particular, the parameter c can be used to modulate
the sparsity of the simulated networks. Each vertex has a label in
$[1,\ldots,q]$, indicating which set it belongs to. A study presented in [29]
identifies conditions under which polynomial-time algorithms can find a
partition that is correlated with the planted partition. More specifically
this is possible when

(6.3) $p_{in}-p_{out}>\frac{\sqrt{qp_{in}+q(q-1)p_{out}}}{\sqrt{n}}.$

This implies that, depending on the generative parameter values, this
condition might not be satisfied. It is often convenient to work with a scaled
affinity matrix so that $p_{in}=\frac{c_{in}}{n}$ and
$p_{out}=\frac{c_{out}}{n}$, where $c_{in}$ is the average vertex degree of
vertices having all edges within the same set and $c_{out}$ is the average
vertex degree of vertices having edges between two different sets. In this
case, the average degree of the graph can be expressed as:

(6.4) $c=\frac{q-1}{q}c_{out}+\frac{1}{q}c_{in}.$

Rewriting ( 6.3) in terms of $c_{in}$ and $c_{out}$, we have:

(6.5) $c_{in}-c_{out}>\sqrt{q[c_{in}+(q-1)c_{out}]},$

and from (6.4) we finally obtain:

(6.6) $c_{in}-c_{out}>q\sqrt{c}.$

We observe that when $c_{in}$ is much greater than $c_{out}$ the graph shows a
strong cluster structure, while such structure becomes weaker when $c_{in}$ is
close to $c_{out}$. If follows that we can implicitly tune the strength of the
set/cluster structure by varying the difference $c_{in}-c_{out}$.

The strength of division of a graph into clusters can be measured by the
modularity [30]. It can be broadly defined as the fraction of the edges that
fall within the given set minus the expected value of such fraction if edges
were distributed at random, in details:

$Q=\frac{1}{2ne}\sum_{ij}[{\mathcal{A}}_{ij}-\frac{c_{i}c_{j}}{2ne}]\delta(C_{i},C_{j}),$

where $ne$ is the number of edges in the whole graph, $c_{i}$, $c_{j}$ are the
degree of vertices $i$ and $j$ respectively,
${\mathcal{A}}=({\mathcal{A}}_{ij})_{ij}$ is the adjacency matrix of the
graph, $\delta(C_{i},C_{j})$ is equal to $1$ if vertices $i$ and $j$ belong to
the same cluster and $0$ otherwise. Graphs with high modularity have dense
connections between vertices within clusters but sparse connections across
clusters. We observe that $Q\in[-1,1]$ and it is often used to evaluate the
quality of a partition obtained by a clustering method. Furthermore, many
clustering methods are based on optimization algorithms which try to maximize
modularity, e.g. Lou and FG.

In some cases it is expected that a graph has natural partition into clusters
referred to as ground truth. Clustering methods aim to estimate such a true
partition. In our simulation framework the ground truth is known and a high
modularity corresponds to $c_{in}>>c_{out}$, while a low modularity
corresponds to $c_{in}\sim c_{out}$. We varied the difference $c_{in}-c_{out}$
above and below the so-called detectability threshold $q\sqrt{c}$, given a
specific set of average degrees $c\in[5,10,15,20]$ and number of clusters
$q\in[2,3,4,5]$. In the following discussion, we use the true labels, i.e.,
the index of the cluster that a vertex belongs to as vertex attribute, so
vertices in the same cluster have the same attribute.

In Figures 4-7, we show a summary of the results obtained for $n=400$
vertices, $q=2,4$, and $c=5,20$. In our simulated data the ground truth labels
are known and we plot the corresponding modularity in subfigure (a) of every
plot. As you can see, our new method $BCMAG$ perfectly matches the truth also
below the detectability threshold, especially for $q=2$. Furthermore, $BCMAG$
largely outperforms $BCMG$, at higher vertex degree and graph density (i.e.,
for $c=20$). On the other hand, it is expected that a method that does not use
the attributes would fail below the detectability threshold. The two methods
$FG$ and $Lou$ do not use labels and maximize the modularity and they both
show an overall overestimation of the modularity in the case of sparse graphs
($c=5$). If we look at the $NMI$, entropy and gain values, we can see that in
all cases for $q=2$, our new method $BCMAG$ performs perfectly well below the
detectability threshold, while for $q=4$, some small discrepancy is observed
when $c_{in}-c_{out}$ is very small (i.e., when modularity is very small).
This analysis shows that $BCMAG$ has very good ability to detect clusters
which are correlated to planted SBM synthetic graphs also when the vertex
groups are not clearly formed and the graph structure is ambiguous. It is
worth noticing that $NMI$ and entropy values for both $Newman$ and $BAGC$, are
not satisfactory especially for sparser graphs ($c=5$) although both methods
make use of the attributes information.

(a)

(b)

(c)

(d)

Figure 4. SBM dataset. Case $n=400$, $q=2$, $c=5$. In each subplot we
indicated the theoretical detectability threshold with a vertical dashed line,
and modularity (a), $NMI$ (b), entropy (c) and gain (d) corresponding to the
different clustering approaches.

(a)

(b)

(c)

(d)

Figure 5. SBM dataset. Case $n=400$, $q=2$, $c=20$. In each subplot we
indicated the theoretical detectability threshold with a vertical dashed line,
and modularity (a), $NMI$ (b), entropy (c) and gain (d) corresponding to the
different clustering approaches.

(a)

(b)

(c)

(d)

Figure 6. SBM dataset. Case $n=400$, $q=4$, $c=5$. In each subplot we
indicated the theoretical detectability threshold with a vertical dashed line,
and modularity (a), $NMI$ (b), entropy (c) and gain (d) corresponding to the
different clustering approaches.

(a)

(b)

(c)

(d)

Figure 7. SBM dataset. Case $n=400$, $q=4$, $c=20$. In each subplot we
indicated the theoretical detectability threshold with a vertical dashed line,
and modularity (a), $NMI$ (b), entropy (c) and gain (d) corresponding to the
different clustering approaches.

###  6.3. LFR-EA graphs

In this section we discuss results obtained on the synthetic graphs introduced
in [31] as reference benchmark to assess performance of clustering algorithms
using vertex attributes information. This benchmark is based on the
Lancichinetti–Fortunato–Radicchi (LFR) model [32], which uses power-law
distributions for both vertex degree and community size. This model generates
more realistic networks with respect to the above SBM, providing heterogeneity
in the distributions of both node degrees and of community sizes. A structure
mixing parameter ranging in $[0,1]$ controls the fraction of edges that are
between communities, i.e., the block structure in the network. If the
structure mixing is equal to $0$, the network has a well partitioned
structure, i.e., it has a large modularity. On the other hand, when the
structure mixing parameter approaches $1$, the structure becomes less modular.
In [31], the authors introduce a model to generate attributes for the LFR
graph vertices. The model parameters are the number $m$ of attributes $l_{j},\
j=1,\ldots,m$, the domain size for the attribute values $n_{j}$ and the
attributes noise level. All the vertices in a cluster are assumed to share the
same attribute domain values. The attributes of some randomly chosen vertices
in each cluster are assigned attribute values affected by the specified noise
level, so that this benchmark model can be used to analyze the robustness of
the algorithms with respect to perturbations in the attribute values. This
benchmark model, named LFR-EA, is often used in the literature [33, 34] to
compare results with those discussed in [31, section 5.2] also because several
methods implementations are not made available by the authors.

We generated LFR-EA with the set of parameters detailed in [31, Table 1]. The
generated graphs have $1000$ vertices and $2$ classes of attributes. Following
the same analysis approach in [31], we considered 9 structure mixing
parameters $[0.1,0.2,\ldots,0.9]$, and 10 attribute noise levels
$[0,0.1,\ldots,0.9]$. For each couple of mixing parameter and attribute noise
level, we generated $100$ samples graphs, for a total of $9000$ test cases. On
each attributed graph, we run a set of clustering methodologies and we
summarised the results in terms of $NMI$ heatmaps in Figure 8. Note that each
color represents the mean $NMI$ value corresponding to $100$ graph samples
generated in correspondence of the specific set of mixing parameter and
attribute noise levels. In each heatmap the x-axis represents the structure
mixing parameter and the y-axis represents the attributes noise.

Figure 8. LFR-EA dataset. Each subplot illustrates the NMI heatmaps for each
method; the color range is based on mean NMI for $100$ graph samples. The
x-axis represents the structure mixing parameter and the y-axis represents the
attributes noise. $BCMG$ (a), $Lou$ (b), $FG$ (c), $BCMAG$ (d) and $BAGC$ (e).

As expected, the clustering algorithms which only use information on the
structure, such as $BCMG$, $Lou$ and $FG$, show the same values of $NMI$ in
vertical columns, indeed attributes have no impact on the performance of those
algorithms. As in the SBM-based graph samples, $BCMG$ confirms his ability to
well identify partitions in highly modular graphs also in the case of LFR
benchmark, and shows better accuracy than $FG$ in recovering the clusters for
decreasing graph modularity. Furthermore, as well known, $Lou$ has the best
performance with a perfect $NMI$ ($\sim 1$) as long as a structure mixing
smaller or equal to $0.5$ is used, and it shows very high values of mean $NMI$
also when the structure mixing reaches the value of $0.6$. However, when
structure mixing parameter increases, all methods are unable to recover
correctly the clusters.

If we look at the average NMI of the methods employing information on
attributes, we observe that increasing attributes noise has more impact on
$BCMAG$ than onto $BACG$, especially when increasing structure mixing
parameter. On the other hand, the top-left and the bottom-right corners of the
$BCMAG$ subplot square show it has larger values of mean $NMI$. This indicates
that the attributes noise have small impact on the ability of $BCMAG$ to
detect the correct clusters when the graph modularity is high (i.e. structure
mixing is low). In the same way, the bottom rows of the $BCMAG$ subplot
indicate that when there is no attributes noise (or the noise level is small),
$BCMAG$ is able to work well also for very small modularities, performing
better than all the other methods. We can conclude that our proposed strategy
$BCMAG$ outperforms $BCMG$ and all the other methods when employing vertices
attributes with no noise, improving the accuracy in recovering the right
clusters also for very large mixing parameter. On the contrary, the impact of
the attributes noise generally increases for small values of mixing parameter.
Note that the plots in [31, Section 5.2, Figure 2 (e)] show that the SA-
cluster method [2, 1], where the graphs attribute augmentation used in this
paper was first introduced, has the worst performance on these benchmark
graphs, not being able to recover the correct clusters. This demonstrates the
benefits of our embedding strategy and of the new vector-valued distance as a
valid alternative to distance-based approach relying on random walks on
attribute-augmented graphs.

###  6.4. Real networks

In this section we show the application of our new method on graphs
originating by two real data networks that we indicate in the following as
$Lazega\ Lawyers$ and $Yeast$. A few indicators for these two networks are
provided in Table 2. In particular we indicate the network names, number of
vertices $|V|$, number of edges $|E|$, sparsity (defined as
$\frac{|E|}{|V|(|V|-1)/2}$), average vertex degree $c$, and number of
attributes. Further details on experimental setting and discussion of results
follow in the subsections.

Name | $|V|$ | $|E|$ | sparsity | $c$ | $m$
---|---|---|---|---|---
Yeast | 2375 | 11693 | 0.004 | 9.85 | 1
adv | 71 | 717 | 0.289 | 20.20 | 7
adv36 | 36 | 289 | 0.459 | 16.06 | 7
friend36 | 36 | 187 | 0.297 | 10.39 | 7
Table 2. Features of real networks. For each network we show the name, number
of vertices $|V|$, number of edges $|E|$, sparsity, average degree $c$ and
number of available attributes $m$

#### Lazega Lawyers dataset

The data $Lazega\ Lawyers$ comes from a network study of corporate law
partnership that was carried out in a Northeastern US corporate law firm
between 1988-1991. It includes measurements of networks among the 71 attorneys
(partners and associates) of the firm. In particular, we looked at the advice
network (adv) and at the friendship network (friend). Various members’
attributes are also part of the dataset, including seniority, formal status,
office in which they work, gender, law school attended, individual performance
measurements (hours worked, fees brought in), attitudes concerning various
management policy options, etc. The ethnography, organizational and network
analyses of this case are available in [35]. The number $36$ indicates that
subset of data for $36$ partners only were considered. In our analysis we
considered 7 attributes, namely: status (1=partner; 2=associate), gender
(1=man; 2=woman), office (1=Boston; 2=Hartford; 3=Providence), firm years,
age, practice (1=litigation; 2=corporate), law school (1: Harvard, Yale; 2:
Ucon; 3: other). The 7 attributes have a different number of distinct values,
as summarized in Table 3.

Attribute | Distinct values
---|---
Status | 2
Gender | 2
Office | 3
Firm years | 27
Age | 33
Practice | 2
Law school | 3
Table 3. Attributes used in the analysis of the $Lazega\ Lawyers$ networks and
the corresponding number of distinct values.

In Figures 9-11, we plot the modularity and the entropy associated to the
methods $BCMAG$, $BCMG$, $BAGC$ and $Lou$, respectively in subfigure (A) and
(B) for different choice of the number $q$ of clusters. The $7$ categories of
attributes considered in the analysis of the 3 $Lazega\ Lawyers$ networks show
a different number of distinct values, and hence there is not an obvious
choice of the number of clusters based on the attributes. As a consequence we
show the number of clusters suggested by the $Lou$ method and the associated
modularity and entropy, computed as in (6.1). As we can see, the modularity of
$BCMAG$ is on average the closest to $Lou$ clustering in every case. On the
other hand, the entropy of $BCMAG$ is the lowest for all cases but the adv36,
where oscillations are observed when $4$ and $6$ clusters are built. In all
cases the new method appears as a large improvement over our previous method
$BCMG$ which does not employ attributes information. It is worth noticing,
that in all cases $BCMAG$ shows the best entropy value when $q=7$, indicating
a good agreement with the partition based on attribute categories.

(a)

(b)

Figure 9. Lazega layers dataset adv. Modularity (a) and entropy (b).

(a)

(b)

Figure 10. Lazega layers dataset adv36. Modularity (a) and entropy (b).

(a)

(b)

Figure 11. Lazega layers dataset friend36. Modularity (a) and entropy (b).

#### Yeast dataset

In this subsection, we apply the described methodologies to a comprehensive
protein-protein interaction network, that we refer to as $Yeast$. This
biological network aims to reveal many aspects of the complex regulatory
network underlying the cellular function. This data set was compiled by [36],
combining various sources. Only the interactions that have “high” and “medium”
confidence are included in the analysis. The data were downloaded from [37]
and are also available in the R package $igraphdata$. We analysed the maximum
connected component using as metadata the single label of the $13$ protein
classes described in [38]. Note that $39$ out of the $2375$ network proteins
had not given class label. We considered each of them as belonging to an
individual class. This resulted in a total of $2375+13+39$ augmented graph
nodes. In Figure 12 we plot the modularity (sub-figure (a)) and the entropy
(sub-figure (b)) corresponding to the clustering provided by $BCMAG$, $BCMG$,
$BAGC$, $Newman$, computed as in (6.2) when ${\mathcal{C}}^{\prime}$ is the
ground truth represented by the attribute labels. We also provide values
obtained by $Lou$ and the attribute labels (as the ground truth), in this case
the entropy is computed as in (6.1). In this example, as said, we have one
single set of labels indicating the protein class. Intuitively, we could
consider as number of clusters $q=13$ (vertical cyan line), corresponding to
the distinct classes provided in the labels. Alternatively, as in the previous
example, we could use the number of clusters $q=23$ found by $Lou$ (vertical
brown line). As a consequence the most interesting area of the plots falls in
the area $13\leq q\leq 23$. As we can see in sub-figure (a), the lowest
modularity corresponds to the protein labels clustering. This indicates that
the attributes provide low separation between the proteins in this case. Note
that the modularity of all the methodologies is less oscillating when $13\leq
q\leq 23$. Moreover, the modularity values provided by $Newman$ and $BCMG$ are
closer to the modularity of $Lou$, while $BCMAG$ and $BAGC$ provide modularity
values that are closer to the labels modularity. In particular, the $BCMAG$
modularity curve offers a trade off between using the labels or not. Indeed
the $BCMAG$ modularity values are on average halfway between the horizontal
cyan line ($Lou$ modularity) and the brown horizontal line (labels
modularity).

If we look at the entropy in sub-figure (b), the general trend of the 4
methodologies $BCMAG$, $BCMG$, $BAGC$ and $Newman$ is to have a decreasing
entropy (conditioned on the true labels) for increasing number of clusters
$q$. The lowest entropy is reached by $BCMAG$ for any $q\in[3,31]$. The
maximum absolute entropy is associated with the labels partition (cyan
horizontal line). As a consequence $BCMAG$ has the maximum gain and hence
appears as the most informative clustering.

(a)

(b)

Figure 12. Yeast dataset. Modularity (a) and entropy (b).

##  7\. Conclusions

In this work we extend a clustering method for undirected graphs to employ
possible additional information, such as vertex attributes, in order to
improve clustering detection and quality. The new method relies on an already
proposed augmentation of the graph which incorporates vertex attributes in the
graph in terms of new vertices and new edges. Starting from this augmentation
we embed the graph in a low-dimensional vector space spanned by some near-null
space vectors computed by an adaptive AMG process approximating small
eigenvalues of the augmented graph Laplacian. Each such vector is
appropriately used to build $1+m$ coordinates, where $m$ is the number of
attributes, of the original graph vertices, so that vertices sharing an
attribute value are in the same hyperplane defined by that value. We also
defined a modified Euclidean distance so that the new method was able to
detect clusters by applying a modified K-means optimization algorithm.
Experimental results show that the new method largely outperforms the previous
approach which ignores the attributes for low-modularity graphs and it is also
better than some of the available methods which fuse attribute and structure
information, also including the SA-cluster which is considered the most
representative of the early-fusion methods. One limitation of our approach can
rely on the fact that it can only be computationally feasible for categorical
attributed networks. Moreover, our approach would require large memory
footprint and computational cost for very large networks, when the number of
distinct attributes values is high. Feature selection or attributes
transformation could be explored to address these issues. Further development
would be needed to handle missing data or sparse attributes.

##  8\. Acknowledgments

The authors thank Clara De Santis, post-graduate fellow at IAC-CNR for helping
in implementing some of the functions dealing with attributed graphs and
included in the code `BCMatch4Graphs`. We also thank Kamal Berahmand and
Haithum Elhadi for sharing with us their codes implementing the LFR-EA model
to generate graph samples.

## References

  * [1] H. Cheng, Y. Zhou, J. Yu, Clustering large attributed graphs: A balance between structural and attribute similarities, ACM Trans. Knowl. Discov. 2 (2011) 12.1–12.33.  doi:10.1145/1921632.1921638.
  * [2] Y. Zhou, H. Cheng, J. Yu, Graph clustering based on structural/attribute similarities, in: Proc. VLDB Endowment, Vol. 2 (1) of VLDB ‘09, August 24-28, 2009, VLDB Endowment, Lyon, France, 2009, pp. 718–729.  doi:10.14778/1687627.1687709.
  * [3] S. Schaeffer, Graph clustering, Computer Science Review 1 (2007) 27–64.  doi:10.1016/j.cosrev.2007.05.001.
  * [4] S. Fortunato, Community detection in graphs, Physics Reports 486 (2010) 75–174.  doi:10.1016/j.physrep.2009.11.002.
  * [5] M. Nascimento, A. de Carvalho, Spectral methods for graph clustering - a survey, European Journal of Operational Research 211 (2011) 221–231.  doi:10.1016/j.ejor.2010.08.012.
  * [6] P. D’Ambra, L. Cutillo, P. Vassilevski, Bootstrap AMG for spectral clustering, Computational and Mathematical Methods 1 (2019) e1020.  doi:10.1002/cmm4.1020.
  * [7] T. Hastie, R. Tibshirani, J. Friedman, The element of statistical learning, Springer, New York, USA, 2001.
  * [8] C. Bothorel, J. D. Cruz, M. Magnani, B. Micenková., Clustering attributed graphs: models, measures and methods, Network Science 3 (3) (2015) 408–444.  doi:10.1017/nws.2015.9.
  * [9] P. Chunaev, Community detection in node-attributed social networks: A survey, Computer Science Review 37, 100286 (2020) 1–24.  doi:10.1016/j.cosrev.2020.100286.
  * [10] P. Chunaev, T. Gradov, B. K., Community detection in node-attributed social networks: how structure-attributes correlation affects clustering quality, Procedia Computer Science 178 (2020) 355–364.  doi:10.1016/j.procs.2020.11.037.
  * [11] Z. Xu, Y. Ke, Y. Wang, H. Cheng, J. Cheng, A model-based approach to attributed graph clustering, in: Proceedings of the 2012 ACM SIGMOD International Conference on Management of Data, SIGMOD ’12, Association for Computing Machinery, New York, NY, USA, 2012, p. 505–516.  doi:10.1145/2213836.2213894.
  * [12] Z. Xu, Y. Ke, Y. Wang, H. Cheng, J. Cheng, Gbagc: A general bayesian framework for attributed graph clustering, ACM Trans. Knowl. Discov. Data 9, 5 (2014) 1–43.  doi:10.1145/2629616.
  * [13] Z. Huang, Y. Ye, X. Li, F. Liu, H. Chen, Joint weighted nonnegative matrix factorization for mining attributed graphs, in: J. Kim, K. Shim, L. Cao, J. Lee, X. Lin, Y. Moon (Eds.), Advances in Knowledge Discovery and Data Mining - 21st Pacific-Asia Conference, PAKDD 2017, Jeju, South Korea, May 23-26, 2017, Proceedings, Part I, Vol. 10234 of Lecture Notes in Computer Science, 2017, pp. 368–380.  doi:10.1007/978-3-319-57454-7\\_29.
URL https://doi.org/10.1007/978-3-319-57454-7_29

  * [14] J. Neville, M. Adler, D. Jensen, Clustering relational data using attribute and link information, in: In Proceedings of the Text Mining and Link Analysis Workshop, 18th International Joint Conference on Artificial Intelligence, 2003, pp. 9–15.
  * [15] D. Combe, C. Largeron, E. Egyed-Zsigmond, M. Géry, Combining relations and text in scientific network clustering, in: First International Workshop on Semantic Social Network Analyis and Design at IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining, 2012, pp. 1280–1285.
  * [16] M. Newman, A. Clauset, Structure and inference in annotated networks, Nature Communications 7:11863 (2016).  doi:10.1038/ncomms11863.
  * [17] U. von Luxburg, A tutorial on spectral clustering, Statistics and Computing 17 (2007) 395–416.  doi:10.1007/s11222-007-9033-z.
  * [18] R. Peng, H. Sun, P. Zanetti, Partitioning well-clustered graphs: spectral clustering works!, Journal of Machine Learning Research: Workshop and Conference Proceedings 40 (2015) 1–33.
  * [19] P. S. Vassilevski, Multilevel block factorization preconditioners, matrix-based analysis and algorithms for solving finite element equations, Springer, New York, USA, 2008.
  * [20] P. D’Ambra, P. Vassilevski, Adaptive AMG with coarsening based on compatible weighted matching., Computing and Visualization in Science 16 (2013) 59–76.  doi:10.1007/s00791-014-0224-9.
  * [21] P. D’Ambra, S. Filippone, P. Vassilevski, BootCMatch: a software package for bootstrap AMG based on graph weighted matching, ACM Trans. Math. Softw. 44-39 (2018) 1–25.  doi:10.1145/3190647.
  * [22] P. D’Ambra, P. Vassilevski, Improving solve time of aggregation-based adaptive AMG, Numerical Linear Algebra with Applications 26 (2019) e2269.  doi:10.1002/nla.2269.
  * [23] V. Kuhlemann, P. Vassilevski, Improving the communication pattern in mat-vec operations for large scale-free graphs by disaggregation, SIAM Journal on Scientific Computing 35 (2013) S465–S486.  doi:10.1137/12088313X.
  * [24] D. V. Blondel, J. L. Guillaume, R. Lambiotte, E. Lefebvre, Fast unfolding of communities in large networks, Journal of Statistical Mechanics Theory and Experiment 2008 (2008) P10008.  doi:10.1088/1742-5468/2008/10/P10008.
  * [25] A. Clauset, M. Newman, C. Moore, Finding community structure in very large networks, Physical Review, E, 70 (2004) 066111.  doi:10.1103/PhysRevE.70.066111.
  * [26] X. Z., Bayesian attributed graph clustering (BAGC code), https://github.com/zhiqiangxu2001/BAGC (2012).
  * [27] L. Danon, J. Duch, A. Diaz-Guilera, A. Arenas, Comparing community structure identification, Journal of Statistical Mechanics: Theory and Experiment 2005 (2005) P09008.  doi:10.1088/1742-5468/2005/09/P09008.
  * [28] P. W. Holland, K. Laskey, S. Leinhardt, Stochastic blockmodels: First steps, Social Networks 5-2 (1983) 109–137.  doi:10.1016/0378-8733(83)90021-7.
  * [29] A. Decelle, F. Krzakala, C. Moor, L. Zdeborová, Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications, Physical Review, E 84 (2011) 066106, 10.1103/PhysRevE.84.066106.  doi:10.1103/PhysRevE.84.066106.
  * [30] M. Newman, Modularity and community structure in networks, Proceedings of the National Academy of Sciences 103 (23) (2006) 8577–8582.  doi:10.1073/pnas.0601602103.
  * [31] H. Elhadi, G. Agam, Structure and attributes community detection: comparative analysis of composite, ensemble and selection methods, in: Proceedings of the 7th Workshop on Social Network Mining and Analysis, SNAKDD ’13, August 2013, ACM, Chicago, Illinois (USA), 2013, pp. 10:1–10:7.  doi:10.1145/2501025.2501034.
  * [32] A. Lancichinetti, S. Fortunato, F. Radicchi, Benchmark graphs for testing community detection algorithms, Phys. Rev. E 78 (2008) 046110.  doi:10.1103/PhysRevE.78.046110.
URL https://link.aps.org/doi/10.1103/PhysRevE.78.046110

  * [33] C. Pizzuti, A. Socievole, Multiobjective optimization and local merge for clustering attributed graphs, IEEE Transactions on Cybernetics 50 (12) (2020) 4997–5009.  doi:10.1109/TCYB.2018.2889413.
  * [34] K. Berahmand, M. Mohammadi, A. Faroughi, R. P. Mohammadiani, A novel method of spectral clustering in attributed networks by constructing parameter-free affinity matrix, Cluster Computing 25 (2022) 869–888.  doi:10.1007/s10586-021-03430-0.
  * [35] E. Lazega, The Collegial Phenomenon: The Social Mechanisms of Cooperation Among Peers in a Corporate Law Partnership, Oxford University Press, Oxford, UK, 2001\.  doi:10.1093/acprof:oso/9780199242726.001.0001.
  * [36] C. von Mering, R. Krause, B. Snel, et al., Comparative assessment of large-scale data sets of protein–protein interactions, Nature 417 (2002) 399–403.  doi:10.1038/nature750.
  * [37] Yeast data set, http://www.nature.com/nature/journal/v417/n6887/suppinfo/nature750.html.
  * [38] Attributes for Yeast data set, https://rdrr.io/cran/igraphdata/man/yeast.html.