# Linear Clustering Process on Networks

Ivan Jokić<EMAIL_ADDRESS>Faculty of Electrical Engineering,
Mathematics, and Computer Science, Delft University of Technology, P.O. Box
5031, 2600 GA Delft, The Netherlands Piet Van Mieghem Faculty of Electrical
Engineering, Mathematics, and Computer Science, Delft University of
Technology, P.O. Box 5031, 2600 GA Delft, The Netherlands

###### Abstract

We propose a linear clustering process on a network consisting of two opposite
forces: attraction and repulsion between adjacent nodes. Each node is mapped
to a position on a one-dimensional line. The attraction and repulsion forces
move the nodal position on the line, depending on how similar or different the
neighbourhoods of two adjacent nodes are. Based on each node position, the
number of clusters in a network, together with each node’s cluster membership,
is estimated. The performance of the proposed linear clustering process is
benchmarked on synthetic networks against widely accepted clustering
algorithms such as modularity, the Louvain method and the non-back tracking
matrix. The proposed linear clustering process outperforms the most popular
modularity-based methods, such as the Louvain method, while possessing a
comparable computational complexity.

##  1 Introduction

Networks [1, 2] abound and increasingly shape our world, ranging from
infrastructural networks (transportation, telecommunication, power-grids,
water, etc.) over social networks to brain and biological networks. In
general, a network consists of a graph or underlying topology and a dynamic
process that takes place on the network. Some examples of processes on a
network are percolation [3] and epidemic spreading [4, 5], that possess a
phase transition [6, 7]. While most real-world processes on networks are non-
linear, linearisation allows for hierarchical structuring of processes on the
network [8].

The identification of communities and the corresponding hierarchical structure
in real-world networks has been an active research topic for decades [9],
although a single, precise definition of a community does not seem to exist
[10, 11]. In network science, a community is defined as a set of nodes that
share links dominantly between themselves, while a minority of links is shared
with other nodes in the network. Newman proposed in [12] a spectral clustering
algorithm that reveals hierarchical structure of a network, by optimising
modularity, a commonly used quality function of a graph partition. A
heuristic, modularity-based two-step clustering algorithm, proposed by Blondel
et al. in [13], has proved to be computationally efficient and performed among
the best in the comparative study conducted in [14]. Recently, Peixoto
proposed in [15] a nested generative model, able to identify nested partitions
at different resolutions, which thus overcomes an existing drawback of a
majority of clustering algorithms, identifying small, but well-distinguished
communities in a large network. Dannon et al. concluded in their comparative
study [16] that those clustering algorithms performing the best tend to be
less computationally efficient. We refer to [9] for a detailed review on
existing clustering algorithms.

Our new idea here is the proposal of a linear clustering process (LCP) on a
graph, where nodes move in one-dimensional space and tend to concentrate in
groups that lead to network communities. Linear means ”proportional to the
graph”, which is needed, because the aim is to cluster the graph and the
process should only help and not distract from our main aim of clustering. A
non-linear process depends intricately on the underlying graph that we want to
cluster and may result in worse clustering! Our LCP leads to a new and non-
trivial graph matrix $W$ in (12) in Theorem 1, whose spectral decomposition is
at least as good as the best clustering result, based on the non-back tracking
matrix [17]. Moreover, the new graph matrix $W$ has a more ”natural” relation
to clustering than the non-back tracking matrix, that was not designed for
clustering initially. Finally, our resulting LCP clustering algorithm seems
surprisingly effective and can compete computationally with any other
clustering algorithm, while achieving generally a better result!

In Section 2, we introduce notations for graph partitioning and briefly review
basic theory on clustering such as modularity and the stochastic block model.
We introduce the linear clustering process (LCP) on network in Section 3,
while the resulting community detection algorithm is described in Section 4
and Section 5. We compare the performance of our LCP algorithm with that of
the non-back tracking matrix, Newman’s and the Louvain algorithm and provide
results in Section 6, after which we conclude.

##  2 Network or graph clustering

A graph $G\left(\mathcal{N},\mathcal{L}\right)$ consists of a set
$\mathcal{N}$ of $N=\left|\mathcal{N}\right|$ nodes and a set $\mathcal{L}$ of
$L=\left|\mathcal{L}\right|$ links and is defined by the $N\times N$ adjacency
matrix $A$, where $a_{ij}=1$ if node $i$ and node $j$ are connected by a link,
otherwise $a_{ij}=0$. The $N\times 1$ degree vector $d$ obeys $d=A\cdot u$,
where the $N\times 1$ all-one vector $u$ is composed of ones. The
corresponding $N\times N$ degree diagonal matrix is denoted by
$\Delta=\text{diag}\left(d\right)$.

The set of neighbours of node $i$ is denoted by $\mathcal{N}_{i}=\\{k\mid
a_{ik}=1,k\in\mathcal{N}\\}$ and the degree of node $i$ equals the cardinality
of that set, $d_{i}=\left|\mathcal{N}_{i}\right|$. The set of common
neighbours of node $i$ and node $j$ is $\mathcal{N}_{i}\cap\mathcal{N}_{j}$,
while the set of neighbours of node $i$ that do not belong to node $j$ is
$\mathcal{N}_{i}\setminus\mathcal{N}_{j}$. The degree of a node $i$ also
equals the sum of the number of common and different neighbours between nodes
$i$ and $j$

$d_{i}=\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|+\left|\mathcal{N}_{i}\cap\mathcal{N}_{j}\right|$
(1)

The number of common neighbours between nodes $i$ and $j$ equals the $ij$-th
element of the squared adjacency matrix

$\left|\mathcal{N}_{i}\cap\mathcal{N}_{j}\right|=\left(A^{2}\right)_{ij}$ (2)

because $(A^{k})_{ij}$ represents the number of walks with $k$ hops between
node $i$ and node $j$ (see [18, p. 32]). From (1), (2) and
$d_{i}=(Au)_{i}=(A^{2})_{ii}$, we have

$\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|=\left(A^{2}\right)_{ii}-\left(A^{2}\right)_{ij}$

and

$\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|+\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|=\left(A^{2}\right)_{ii}+\left(A^{2}\right)_{jj}-2\left(A^{2}\right)_{ij}$

The latter expression is analogous to the effective resistance $\omega_{ij}$
between node $i$ and node $j$,

$\omega_{ij}=Q_{ii}^{{\dagger}}+Q_{jj}^{{\dagger}}-2Q_{ij}^{{\dagger}}$

in terms of the pseudoinverse $Q_{ii}^{{\dagger}}$ of the Laplacian matrix
$Q=$ $\Delta-A$ (see e.g. [19]).

Network clustering aims to organize nodes in the graph $G$ in such a way that
the majority of links in the graph connect nodes of the same cluster, while a
minority of links is shared by nodes from different clusters. This general
description is imprecise and part of the complication, which causes that
community detection in a network is a challenging problem, that has been
actively researched over the last few decades [9].

Before introducing our linear clustering process (LCP) in Section 3, we
briefly present basic graph partitioning concepts, while the overview of the
more popular clustering methods is deferred to Appendix A.

###  2.1 Network modularity

Newman and Girvan [20] proposed the modularity as a concept for a network
partitioning,

$m=\frac{1}{2L}\cdot\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}\left(a_{ij}-\frac{d_{i}\cdot
d_{j}}{2L}\right)\cdot\mathbf{1}_{\\{i\,\text{and}\,j\,\in\text{ same
cluster}\\}},$ (3)

where $\mathbf{1}_{x}$ is the indicator function that equals 1 if statement
$x$ is true, otherwise $\mathbf{1}_{x}=0$. The modularity $m$ compares the
number of links between nodes from the same community with the expected number
of intra-community links in a network with randomly connected nodes. When the
modularity $m$ close to $0$, the estimated partition is as good as a random
partition would be. On the contrary, a modularity $m$ close to $1$ indicates
that the network can be clearly partitioned into clusters. Defining the
$N\times N$ modularity matrix $C$,

$C_{ij}=\begin{cases}1&\text{if nodes $i$ and $j$ belong to the same
cluster}\\\ 0&\text{otherwise},\end{cases}$ (4)

allows us to rewrite the modularity (3) as a quadratic form,

$m=\frac{1}{2L}\cdot u^{T}\cdot\left(A\circ C-\frac{1}{2L}\cdot\left(d\cdot
d^{T}\right)\circ C\right)\cdot u,$ (5)

where $\circ$ denotes the Hadamard product [21]. The number of clusters in a
network is denoted by $c$, while the $c\times 1$ vector $n$ defines the size
of each cluster:

$n=\begin{bmatrix}n_{1}&n_{2}&\dots&n_{c}\end{bmatrix},$ (6)

where the number of nodes in cluster $i$ is denoted as $n_{i}$.

###  2.2 Stochastic block model

The performance of the clustering methods in this paper are benchmarked on
random graphs, generated by the Stochastic Block Model (SBM), proposed by
Holland [22]. The SBM model generates a random graph with community structure,
where a link between two nodes exists with different probability, depending on
whether the nodes belong to the same cluster or not.

In this paper, we focus on the symmetric stochastic block model (SSBM), where
only two different link probabilities are defined. Two nodes are connected via
a link with probability $p_{in}$ if they belong to the same cluster, otherwise
the direct link exists with probability $p_{out}$. Communities emerge when the
link density within clusters is larger than the inter-community link
probability $p_{in}>p_{out}$. Furthermore, we restrict clusters to be of the
same size:

$n_{i}=\frac{N}{c},i\in\\{1,\,2,\,\dots\,,c\\},$

which causes the expected degree to be the same for each node:

$E[D]=\frac{b_{in}+\left(c-1\right)\cdot b_{out}}{c},$ (7)

irrespective of its cluster membership. We further consider a sparse and
assortative SSBM. The SBMM is called sparse and assortative when the link
probabilities $p_{in}=\frac{b_{in}}{N}$ and $p_{out}=\frac{b_{out}}{N}$ are
defined upon positive constants $b_{in}>b_{out}$ that stay constant when
$N\rightarrow\infty$. Decelle et al. [23, 24] found that when the difference
$b_{in}-b_{out}$ is above the detectability threshold

$b_{in}-b_{out}>c\cdot\sqrt{E[D]},$ (8)

it is theoretically possible to recover cluster membership of the nodes,
otherwise, the community structure of a network is not distinguishable from
randomness. The threshold (8) marks a phase transition between the
undetectable and the theoretically detectable regime of the SSBM.

##  3 Linear clustering process (LCP) on a graph

###  3.1 Concept of the clustering process

Each node $i$ in the graph $G$ is assigned a position $x_{i}[k]$ on a line
(i.e. in one-dimensional space) at discrete time $k$. We define the $N\times
1$ position vector $x[k]$ at discrete time $k$, where the $i$-th vector
component consists of the position $x_{i}[k]$ of node $i$ at time $k$. We
initialize the $N\times 1$ position vector $x[0]$ by placing nodes
equidistantly on the line and assign integer values from $1$ to $N$ to the
nodes, thus, $x[0]=\begin{bmatrix}1&2&\dots&N\end{bmatrix}^{T}$. At last, we
restrict the position $x_{i}[k]$ to positive real values.

We propose a dynamic process that determines the position of nodes over time.
The position difference between nodes of the same cluster is relatively small.
On the contrary, nodes from different clusters are relatively far away, i.e.
their position difference is relatively high. Based on the position vector
$x[k]$, we will distinguish clusters, also called communities, in the graph
$G$.

The proposed clustering process consists of two opposite and simultaneous
forces that change the position of nodes at discrete time $k$:

Attraction

Adjacent nodes sharing many neighbours are mutually attracted with a force
proportional to the number of common neighbours. In particular, the attractive
force between node $i$ and its neighboring node $j$ is proportional to
$\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)$,
where $\alpha$ is the attraction strength and
$\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)$ equals the
number of common neighbors plus the direct link, i.e. $a_{ij}=1$.

Repulsion

Adjacent nodes sharing a few neighbours are repulsed with a force proportional
to the number of different neighbours. The repulsive force between node $i$
and its neighboring node $j$ is proportional to
$\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|-1\right)$,
where $\delta$ is the repulsive strength and
$\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|-1\right)$ equals
the set of neighbours of node $j$ that do not belong to node $i$ minus the
direct link (that is included in
$\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|)$. Since the force
should be symmetric and the same if $i$ and $j$ are interchanged, we end up
with a resultant repulsive force proportional to
$\frac{1}{2}\cdot\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)$.

###  3.2 LCP in discrete time

Since computers operate with integers and truncated real numbers, we
concentrate on discrete-time modeling. The continuous-time description is
derived in Appendix B. We denote the continuous-time variables by
$y\left(t\right)$ and the continuous time by $t$, while the discrete-time
counterpart is denoted by $y\left[k\right]$, where the integer $k$ denotes the
discrete time or $k$-th timeslot. The transition from the continuous-time
derivative to the discrete-time difference is

$\displaystyle\frac{dx_{i}(t)}{dt}=\lim_{\Delta t\rightarrow
0}\frac{x_{i}(t+\Delta t)-x_{i}(t)}{\Delta t}\rightarrow$
$\displaystyle\left.\frac{x_{i}(t+\Delta t)-x_{i}(t)}{\Delta t}\right|_{\Delta
t=1}\overset{\text{def}}{=}x_{i}\left[k+1\right]-x_{i}\left[k\right]$

Corresponding to the continuous-time law in Appendix 38 and choosing the time
step $\Delta t=1$, the governing equation of position $x_{i}[k]$ of node $i$
at discrete time $k$ is

$x_{i}[k+1]=x_{i}[k]+\sum\limits_{j\in\mathcal{N}_{i}}\bigg{(}\frac{\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}{d_{j}d_{i}}-\frac{\frac{1}{2}\cdot\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{d_{j}d_{i}}\bigg{)}\cdot\Big{(}x_{j}[k]-x_{i}[k]\Big{)}$
(9)

where $\alpha$ and $\delta$ are, in the discrete-time setting, the strength
(in dimensionless units) for attraction and repulsion, respectively. The
maximum position difference at the initial state is $x_{N}[0]-x_{1}[0]=N-1$.

Node $j$ attracts an adjacent node $i$ with force proportional to their
position difference ($x_{j}[k]-x_{i}[k]$). The intensity of the attractive
force decreases as nodes $i$ and $j$ are closer on a line. The attraction is
also proportional to the number common neighbours
$\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|$ of node $i$ and node $j$
plus the direct link, as nodes tend to share most links with other nodes from
the same cluster. On the contrary, node $j$ repulses node $i$ with a rate
proportional to their position difference ($x_{j}[k]-x_{i}[k]$) and the
average of the number of node $j$ neighbours
$\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|$ that are not connected
to the node $i$ and, similarly, the number of node $i$ neighbors,
$\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|$ that are not connected
to the node $j$. The repulsive and attractive force are, as mentioned above,
symmetric in strength, but opposite, if $i$ is interchanged by $j$.

Figure 1: Dependence of the attractive and repulsive force on the number of
common neighbours of adjacent nodes $i$ and $j$ (left-figure). Directions of
the attraction and repulsion forces between the adjacent nodes (middle-
figure). Dependence of the attractive and repulsive force on the absolute
position distance between adjacent nodes $i$ and $j$ (right-figure).

The directions of both attractive and repulsive forces between two adjacent
nodes $i$ and $j$ as well the dependence of both forces on the number of
common neighbours $\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|$ and the
absolute position distance $\left|x_{j}[k]-x_{i}[k]\right|$ are illustrated in
Fig 1.

In the continuous-time setting, as provided in Appendix 38, we eliminate one
parameter by scaling the time $t^{\ast}=\delta t$. Because the time step
$\Delta t=1$ is fixed and cannot be scaled, the discrete-time model consists
of two parameters $\alpha\geq 0$ and $\delta\geq 0$.

So far, we have presented an additive law, derived in the common Newtonian
approach. The corresponding multiplicative law in discrete time is

$x_{i}[k+1]=x_{i}[k]\cdot\Bigg{(}1+\sum\limits_{j\in\mathcal{N}_{i}}\bigg{(}\frac{\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}{d_{i}\cdot
d_{j}}-\frac{\frac{1}{2}\cdot\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{d_{i}\cdot
d_{j}}\bigg{)}\cdot\Big{(}x_{j}[k]-x_{i}[k]\Big{)}\Bigg{)}$

(10)

Although the physical intuition is similar, the multiplicative process in (10)
behaves different in discrete time than the additive law in (9). Since also
the analysis is more complicated, we omit a further study of the
multiplicative law.

We present the analogon of (9) in a matrix form:

###### Theorem 1

The discrete time process (9) satisfies the linear matrix difference equation

$x[k+1]=\left(I+W-\text{diag}\left(W\cdot u\right)\right)\cdot x[k],$ (11)

where the $N\times 1$ vector $u$ is composed of ones, the $N\times N$ identity
matrix is denoted by $I$, while the $N\times N$ topology-based matrix $W$ is
defined as

$W=\left(\alpha+\delta\right)\Delta^{-1}\cdot\left(A\circ
A^{2}+A\right)\cdot\Delta^{-1}-\frac{1}{2}\cdot\delta\left(\Delta^{-1}\cdot
A+A\cdot\Delta^{-1}\right)$ (12)

where $\circ$ denotes the Hadamard product. In particular,

$w_{ij}=a_{ij}\frac{\alpha\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)-\delta\left(\frac{\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|}{2}-1\right)}{d_{i}d_{j}}$
(13)

The explicit solution of the difference equation (11) is

$x[k]=\left(I+W-\text{diag}\left(W\cdot u\right)\right)^{k}x[0]$ (14)

where the $k$-th component of the initial position vector is $(x[0])_{k}=k$.

Proof: Appendix C.1.

Theorem 1 determines the position of the nodal vector $x[k]$ at time $k$ and
shows convergence towards a state, where the sum of attractive and repulsive
forces (i.e. the resulting force) acting on a node are in balance. Nodes with
similar neighbourhoods are grouped on the line, i.e. in the one-dimensional
space, while nodes with a relatively small number of common neighbours are
relatively far away. A possible variant of the proposed linear clustering
process may map the nodal position into a higher dimensional space, like a
circular disk or square in two dimensions, and even with a non-Euclidean
distance metric.

###  3.3 Time-dependence of the linear clustering process

The $N\times N$ matrix $I+W-\text{diag}\left(W\cdot u\right)$ in the governing
equation (11) has interesting properties. As shown in this section, the
related matrix $W-\text{diag}\left(W\cdot u\right)$ belongs to the class of
$M$-matrices, whose eigenvalues have a non-negative real part. The (weighted)
Laplacian is another element of the $M$-matrix class.

###### Property 1

The matrix $I+W-\text{diag}\left(W\cdot u\right)$ is a non-negative matrix.

Proof: The governing equation (11)

$x[k+1]=\left(I+W-\text{diag}\left(W\cdot u\right)\right)\cdot x[k]$

holds for any non-negative vector $x[k]$. Let $x[0]=e_{m}$, the basic vector
with components $\left(e_{m}\right)_{i}=\delta_{mi}$ and $\delta_{mi}$ is the
Kronecker delta, then we find that the $m$-th column

$x[1]=\left(I+W-\text{diag}\left(W\cdot u\right)\right)_{col(m)}$

must be a non-negative vector. Since we can choose $m$ arbitrary, we have
established that $I+W-\text{diag}\left(W\cdot u\right)$ is a non-negative
matrix. $\hfill\square$

###### Property 2

The principal eigenvector of the matrix $I+W-\text{diag}\left(W\cdot u\right)$
is the all-one vector $u$ belonging to eigenvalue 1. All other eigenvalues of
matrix $I+W-\text{diag}\left(W\cdot u\right)$ are real and, in absolute value,
smaller than 1.

Proof: Appendix C.2.

The linear discrete-time system in (11) converges to a steady-state, provided
that
$\lim_{k\rightarrow\infty}||x[k+1]||=\lim_{k\rightarrow\infty}||x[k]||=||x_{s}||$,
which is only possible if the matrix $\left(I+W-\text{diag}\left(W\cdot
u\right)\right)$ has all eigenvalues in absolute value smaller than 1 and the
largest eigenvalue is precisely equal to 1\. Property 2 confirms convergence
and indicates that the steady-state vector $x_{s}=u$ in which the position of
each node is the same. However, the steady state solution $x_{s}=u$ is a
trivial solution, as observed from the governing equation in (9), because the
sum vanishes and the definition of the steady state tells that $x[k+1]=x[k]$,
which is obeyed by any discrete-time independent vector. In other words, the
matrix equation (11) can be written as

$x[k+1]-x[k]=\left(W-\text{diag}\left(W\cdot
u\right)\right)\cdot\left(x[k]-u\right)$

which illustrates that, if $x[k]$ obeys the solution, then $r[k]=x[k]+s\cdot
u$ for any complex number $s$ is a solution, implying that a shift in the
coordinate system of the positions does not alter the physics.

Let us denote the eigenvector $y_{k}$ belonging to the $k$-th eigenvalue
$\beta_{k}$ of the matrix $W-\text{diag}\left(W\cdot u\right)$, where
$\beta_{1}\geq\beta_{2}\geq\cdots\geq\beta_{N}$, then the eigenvalue
decomposition of the real, symmetric matrix is

$W-\text{diag}\left(W\cdot u\right)=Y\text{diag}(\beta)Y^{T}$

where the eigenvalue vector $\beta=(\beta_{1},\beta_{2},\cdots,\beta_{N})$ and
$Y$ is the $N\times N$ orthogonal matrix with the eigenvectors
$y_{1},y_{2},\cdots,y_{N}$ in the columns obeying $Y^{T}Y=YY^{T}=I$. Since
$\beta_{1}=0$ and $y_{1}=\frac{u}{\sqrt{N}}$, it holds for $k>1$ that
$u^{T}y_{k}=0$, which implies that the sum of the components of eigenvector
$y_{k}$ for $k>1$ is zero (just as for any weighted Laplacian [19]). The
position vector in (14) is rewritten as

$x[k]=Y\text{diag}(1+\beta)^{k}Y^{T}x[0]=\sum_{j=1}^{N}(1+\beta_{j})^{k}y_{j}\left(y_{j}^{T}x[0]\right)$

Hence, we arrive at

$x[k]-\frac{u^{T}x[0]}{\sqrt{N}}u=\sum_{j=2}^{N}(1+\beta_{j})^{k}(y_{j}^{T}x[0])\;y_{j}$
(15)

As explained above, the left-hand side is a translated position vector and
physically not decisive for the clustering process. Since $-1<\beta_{j}<0$ for
$j>1$, relation (15) indicates that, for $k\rightarrow\infty$, the right-hand
side tends to zero and the steady-state solution is clearly uninteresting for
the clustering process. We rewrite (15) as

$x[k]-\frac{u^{T}x[0]}{\sqrt{N}}u=(1+\beta_{2})^{k}\bigg{(}(y_{2}^{T}x[0])\;y_{2}+\sum_{j=3}^{N}\left(\frac{1+\beta_{j}}{1+\beta_{2}}\right)^{k}(y_{j}^{T}x[0])\;y_{j}\bigg{)}.$

Since $|1+\beta_{2}|>|1+\beta_{3}|$, we observe that

$\frac{x[k]-\frac{u^{T}x[0]}{\sqrt{N}}u}{(1+\beta_{2})^{k}\;(y_{2}^{T}x[0])}=y_{2}+O\left(\frac{1+\beta_{3}}{1+\beta_{2}}\right)^{k},$
(16)

which tells us that the left-hand side, which is a normalized or scaled,
shifted position vector, tends to the second eigenvector $y_{2}$ with an error
that exponentially decreases in $k$. Hence, for large enough $k$, but not too
large $k$, the scaled shifted position vector provides us the information on
which we will cluster the graph.

The steady state in Property 2 can be regarded as a reference position of the
nodes and does not affect the LCP process nor the $N\times 1$ eigenvector
$y_{2}$, belonging to the second largest eigenvalue $(1+\beta_{2})$ of the
$N\times N$ “operator” matrix $I+W-\text{diag}(W\cdot u)$, which is analogous
to Fiedler clustering based on the $N\times N$ Laplacian $Q$. While the
Laplacian matrix $Q$ essentially describes diffusion and not clustering, our
operator $I+W-\text{diag}(W\cdot u)$ changes the nodal positions, based on
attraction and repulsion, from which clustering naturally arises.

###  3.4 Influence of $\alpha$ and $\delta$ on the eigenvalues $\beta_{k}$ and
the eigenvector $y_{2}$

###### Property 3

The two parameters in the matrix $W$ in (12) satisfy the bounds

$\displaystyle 0\leq\alpha$ $\displaystyle\leq$
$\displaystyle\frac{d_{\max}-1}{d_{\max}-\frac{1}{2}\left(1+\frac{d_{\min}}{d_{\max}}\right)}\leq
1$ (17) $\displaystyle 0\leq\delta$ $\displaystyle\leq$
$\displaystyle\frac{1}{d_{\max}-\frac{1}{2}\left(1+\frac{d_{\min}}{d_{\max}}\right)}$
(18)

Proof: Appendix C.3.

Figure 2: Gap $\beta_{2}-\beta_{3}$ between the second and the third largest
eigenvalue of the $N\times N$ matrix $W-\text{diag}\left(W\cdot u\right)$, for
different values of the attractive $\alpha$ and repulsive $\delta$ strength
(left figure). Relative difference $\frac{\beta_{3}-\beta_{2}}{\beta_{3}}$,
for different values of the attractive $\alpha$ and repulsive $\delta$
strength (right figure). An SSBM network with $N=1000$ nodes, $c=5$ clusters
is used for both plots, with $b_{in}=25$ and $b_{out}=2.5$.

Figure 2 shows that influence of the attractive and repulsive strength
$\alpha$ and $\delta$ on the eigenvalue gap $\beta_{2}-\beta_{3}$ is
relatively small if $\alpha$ and $\delta$ are not too small and obeying the
bounds (17) and (18). While the difference increases when the attraction
strength $\alpha$ is increasing, the repulsive strength $\delta$ has no
visible influence on the eigenvalue gap.

The eigenvalue $\beta_{2}$ depends on the community structure of a graph.
Figure 3 reveals positive correlation between the eigenvalue $\beta_{2}$ and
the modularity index $m$ of a graph. As the modularity index increases, the
eigenvalue $\beta_{2}$ approaches value $1$. In the limit case, when there are
only intra-community links in the network, $\beta_{2}=1$, indicating the
eigenvector $y_{2}$ represents a steady state.

Figure 3: The eigenvalue $\beta_{2}$ versus the modularity index $m$ of an
SSBM graph of $N=999$ nodes and $c=3$ clusters, and an SSBM graph of $N=1000$
nodes, with $c=4,10$ clusters, respectively. The parameters $b_{in}$ and
$b_{out}$ are varied, while keeping average degree $d_{av}=7$ fixed. For each
combination of $b_{in}$ and $b_{out}$, the modularity index $m$ and the
eigenvalue $\beta_{2}$ are computed. The correlation is presented in case only
interactions between direct neighbours are allowed (left figure) and in case
interactions between each pair of nodes are allowed (right figure).

Figure 4 reveals that the repulsive strength $\delta$ does not affect the
eigenvector $y_{2}$ components significantly. Eigenvector $y_{2}$ components
of nodes from the same cluster are better distinguished from the remaining
components of $y_{2}$ for smaller values of repulsive strength $\delta$.

Figure 4:  Sorted eigenvector $\hat{y_{2}}$ components for different values of
the repulsive strength $\delta$, in case of a SSBM network of $N=100$ nodes,
$c=4$ clusters and with parameters $b_{in}=25$ and $b_{out}=1$. The attraction
rate equals $\alpha=0.05$ (left figure) and $\alpha=0.5$ (right figure), while
the repulsive strength $\delta$ obeys bounds in (18).

##  4 From the eigenvector $y_{2}$ to clusters in the network

The interplay of the attractive and repulsive force between nodes drives the
nodal position in discrete time $k$ eventually towards a steady state
$\lim\limits_{k\to\infty}x[k]=u$. However, the scaled and shifted position
vector $x[k]$ in (16) converges in time towards the second eigenvector $y_{2}$
with an exponentially decreasing error. In this section, we estimate the
clusters in network, based on the eigenvector $y_{2}$.

By sorting the eigenvector $y_{2}$ to $\hat{y_{2}}$, the components of $y_{2}$
are reordered and the corresponding relabeling of the nodes of the network
reveals a block diagonal structure of the adjacency matrix $A$. We define the
$N\times N$ permutation matrix $R$ in a way the following equalities hold:

$\displaystyle\hat{y_{2}}$ $\displaystyle=R\cdot y_{2},$ (19)
$\displaystyle\hat{(y_{2})_{i}}=\left(y_{2}\right)_{r_{i}}$
$\displaystyle\leq\hat{(y_{2})_{j}}=\left(y_{2}\right)_{r_{j}},$
$\displaystyle i<j,$

where the $N\times 1$ ranking vector $r=R\cdot w$ and $w=[1,2,\dots,N]$, with
$r_{i}$ denoting the node $i$ ranking in the eigenvector $y_{2}$. The
permutation matrix $R$ allow us to define the $N\times N$ relabeled adjacency
matrix $\hat{A}$, the $N\times 1$ relabeled degree vector $\hat{d}$ of $G$,
and the $N\times 1$ sorted eigenvector $\hat{y_{2}}$ as follows:

$\begin{cases}\hat{A}&=R^{T}\cdot A\cdot R\\\ \hat{d}&=R\cdot d\\\
\hat{y_{2}}&=R\cdot y_{2}.\end{cases}$ (20)

Figure 5: Adjacency matrix $A$ of an SSBM network of $N=1000$ nodes, $c=4$
clusters and parameters $b_{in}=26$, $b_{out}=0.67$ (top-left). Eigenvector
$y_{2}$ components (top-right). Sorted eigenvector $\hat{y_{2}}$ components
(bottom-right). Relabeled adjacency matrix $\hat{A}$ based on the sorted
eigenvector $\hat{y_{2}}$ (bottom-left).

Groups of nodes that have relatively small difference in the eigenvector
$y_{2}$ components, while relatively large difference compared to other nodes
in the network, compose a cluster. Therefore, the community detection problem
transforms into recognizing intervals of similar values in the sorted
eigenvector $\hat{y_{2}}$.

Figure 5 exemplifies the idea, where the adjacency matrix $A$ of a randomly
labeled SSBM network of $N=1000$ nodes and $c=4$ clusters is presented in the
upper-left part, as a heat map. The eigenvector $y_{2}$ is drawn in the upper-
right part, while the sorted eigenvector $\hat{y_{2}}$ is drawn on the bottom-
right side. Finally, the relabeled adjacency matrix $\hat{A}$, based on nodal
ranking of $y_{2}$ is depicted on the lower-left side. The sorted eigenvector
$\hat{y_{2}}$ reveals a stair with four segments, equivalent to four block
matrices on the main diagonal in relabeled adjacency matrix $\hat{A}$.

The eigenvector $y_{2}$ represents a continuous measure of how similar
neighbours of two nodes are. There are two different approaches to identify
network communities for a given eigenvector $y_{2}$:

  * •

Cluster identification based on the sorted eigenvector $\hat{y_{2}}$. This
approach is explained in subsection 4.1.

  * •

Cluster identification based on the ranking vector $r$. This approach does not
rely on the eigenvector $y_{2}$ components, but solely on nodal ranking, as
explained in subsection 4.2.

###  4.1 Community detection based on nodal components of the eigenvector
$y_{2}$

To identify clusters, we observe the difference in eigenvector $y_{2}$
components between nodes with adjacent ranking. If
$\left(\hat{y_{2}}\right)_{i+1}-\left(\hat{y_{2}}\right)_{i}<\theta$, where
$\theta$ denotes a predefined threshold, then the nodes $r_{i}$ and $r_{i+1}$
belong to the same cluster, else the nodes $r_{i}$ and $r_{i+1}$ are
boundaries of two adjacent clusters. The resulting cluster membership function
is

$C_{r_{i+1},r_{i}}=\begin{cases}1&\text{$\left(\hat{y_{2}}\right)_{i+1}-\left(\hat{y_{2}}\right)_{i}<\theta$}\\\
0&\text{otherwise},\end{cases}$ (21)

where the threshold value $\theta$ is determined heuristically, as provided in
section 6. The cluster estimation in (21) can be improved by using other more
advanced approaches, such as the K-means algorithm.

###  4.2 Modularity-based community detection

By implementing (4) and (20) into (3) we obtain:

$m=\frac{1}{2L}\cdot
u^{T}\cdot\left(\hat{A}\circ\hat{C}-\frac{1}{2L}\cdot\left(\hat{d}\cdot\hat{d}^{T}\right)\circ\hat{C}\right)\cdot
u,$ (22)

where $\hat{C}=R^{T}\cdot C\cdot R$. As shown in Figure 5, the network
relabeling based on the ranking vector $r$ reveals block diagonal structure in
$\hat{A}$. Thus, the relabeled modularity matrix $\hat{C}$ has the following
block diagonal structure:

$\hat{C}=\begin{bmatrix}J_{n_{1}\times n_{1}}&O_{n_{1}\times
n_{2}}&\dots&O_{n_{1}\times n_{c}}\\\ O_{n_{2}\times n_{1}}&J_{n_{2}\times
n_{2}}&\dots&O_{n_{1}\times n_{c}}\\\ \vdots&\vdots&\dots&\vdots\\\
O_{n_{1}\times n_{1}}&O_{n_{c}\times n_{2}}&\dots&J_{n_{c}\times n_{c}}\\\
\end{bmatrix},$ (23)

where $c$ denotes number of clusters in network, where the $i$-th cluster is
composed of $n_{i}$ nodes. We define the $N\times 1$ vectors $\hat{e_{i}}$ for
$i=\\{1,2,\dots,c\\}$ as

$\hat{e_{i}}=\begin{bmatrix}O_{\left(1\times\sum_{j=1}^{i-1}n_{j}\right)}&u_{\left(1\times
n_{i}\right)}&O_{\left(1\times\sum_{j=i+1}^{N}n_{j}\right)}\end{bmatrix}^{T},$
(24)

that allows us to redefine
$\hat{C}=\sum_{i=1}^{c}\hat{e_{i}}\cdot\hat{e_{i}}^{T}$ and further simplify
(22):

$m=\frac{1}{2L}\cdot\sum\limits_{i=1}^{c}\hat{e_{i}}^{T}\cdot\left(\hat{A}-\frac{1}{2L}\cdot\left(\hat{d}\cdot\hat{d}^{T}\right)\right)\cdot\hat{e_{i}}.$
(25)

Since the vector $\hat{e_{i}}$ consists of zeros and ones, the equation (25)
represents the sum of elements of the matrix
$\left(\hat{A}-\frac{1}{2L}\cdot\left(\hat{d}\cdot\hat{d}^{T}\right)\right)$
corresponding to each individual cluster.

Algorithm 1 Recursive algorithm for cluster estimation

1:$\hat{A}$ and $\hat{d}$ are the relabeled adjacency matrix $A$ and the
degree vector $d$ (20), while $L$ denotes number of links. The modularity
threshold is denoted by $\theta$. The function returns the $c\times 1$ vector
$b$, whose elements are cluster borders in a relabeled graph.

2:

3:function EstimateClusters($\hat{A},\,\hat{d},\,N,\,L,\,\theta$)

4: $d_{f},d_{b},p,q\leftarrow O_{N\times 1}$

5: $\left(d_{f}\right)_{1}\leftarrow\hat{d}_{1}$

6: $\left(d_{b}\right)_{N}\leftarrow\hat{d}_{N}$

7: $p_{1}\leftarrow-\frac{\hat{d}_{1}^{2}}{(2L)^{2}}$

8: $q_{N}\leftarrow-\frac{\hat{d}_{N}^{2}}{(2L)^{2}}$

9: for $i\leftarrow 2\textrm{ to }N$ do

10: $l\leftarrow N-i$

11: $\left(d_{f}\right)_{i}\leftarrow\left(d_{f}\right)_{i-1}+\hat{d}_{i}$

12:
$\left(d_{b}\right)_{N-i+1}\leftarrow\left(d_{b}\right)_{N-i+2}+\hat{d}_{N-i+1}$

13:
$s\leftarrow\frac{\sum_{j=1}^{i}\hat{a}_{ij}}{L}-\frac{2\cdot\hat{d}_{i}\cdot(d_{f})_{i-1}+\hat{d}_{i}^{2}}{(2L)^{2}}$

14:
$t\leftarrow\frac{\sum_{j=1}^{i}\hat{a}_{N-j+1,l+1}}{L}-\frac{2\cdot\hat{d}_{l+1}\cdot(d_{b})_{l+2}+\hat{d}_{l+1}^{2}}{(2L)^{2}}$

15: $p_{i}\leftarrow p_{i-1}+s$

16: $q_{l+1}\leftarrow q_{l+2}+t$

17: end for

18: $r\leftarrow\operatorname*{arg\,max}_{\mathcal{N}}\left(p+q\right)$

19: if $(p+q)_{r}>\theta$ then

20: $\hat{A_{1}},\,\hat{d_{1}},\,N_{1}\leftarrow$ sub-matrix(vector)
corresponding to the first cluster $\\{1,2,\dots,r\\}$

21: $\hat{A_{2}},\,\hat{d_{2}},\,N_{2}\leftarrow$ sub-matrix(vector)
corresponding to the second cluster $\\{r+1,r+2,\dots,N\\}$

22: return
$\hat{b}\leftarrow\begin{bmatrix}\text{EstimateClusters(}\hat{A_{1}},\hat{d_{1}},N_{1},L,p_{r}\text{)}\\\
r,\\\
\text{EstimateClusters(}\hat{A_{2}},\hat{d_{2}},N_{2},L,q_{r}\text{)}\end{bmatrix}$

23: else

24: return $\hat{b}\leftarrow\emptyset$

25: end if

26:end function

We estimate clusters for a given ranking vector $r$ by optimising the
modularity $m$ recursively. In the first iteration, we examine all possible
partitions of the network in two clusters and compute their modularity. The
partition that generates the highest modularity is chosen. In the second
iteration, we repeat for each subgraph the same procedure and find the best
partitions into two clusters. Once we determine the best partitions for both
subgraphs, we adopt them if the obtained modularity of the generated partition
exceeds the modularity of a parent cluster from the previous iteration. The
recursive procedure stops when the modularity $m$ cannot be further improved,
as described by pseudocode (1). This version of the proposed process is
denoted as LCP in section 6.

###  4.3 Modularity-based community detection for a known number of
communities

The algorithm 1 also applies for graph partition with a known number of
communities $c$. In that case, instead of stopping the recursive procedure
described in algorithm 1 when the modularity $m$ cannot be further improved,
we stop at iteration $(\log_{2}c+1)$. In each iteration, the partition with
the maximum modularity is accepted, even if negative.

As a result, we obtain $2c$ estimated clusters with the $2c\times 2c$
aggregated modularity matrix $M_{c}$:

$\left(M_{c}\right)_{gh}=\sum\limits_{i\in g,j\in
h}\left(\hat{A}-\frac{1}{2L}\cdot\hat{d}\cdot\hat{d}^{T}\right)_{ij},$ (26)

where $g,h\in\\{1,2,\dots,2c\\}$ denote estimated communities. The aggregated
modularity matrix $M_{c}$ allows us to merge adjacent clusters, until we reach
$c$ communities in an iterative way. We observe the $(2c-1\times 1)$ vector
$\mu$, where $\mu_{g}=\left(M_{c}\right)_{g,g+1}$. The maximum element of
$\mu$ indicates which two adjacent clusters can be merged, so that modularity
index $m$ is negatively affected the least. By repeating this procedure $c$
times, we end up with the graph partition in $c$ clusters. This version of the
proposed process is denoted as LCPn in Section 6.

###  4.4 Non-back tracking method versus LCP

Angel et al. [25, p.12] noted that the $2N$ non-trivial eigenvalues of the
$2L\times 2L$ non-back tracking matrix $B$ from (37) are contained in
eigenvalues of the $2N\times 2N$ matrix $B^{*}$:

$B^{*}=\begin{bmatrix}A&I-\Delta\\\ I&O\end{bmatrix},$ (27)

where the $N\times N$ matrix with all zeros is denoted as $O$. The $2N\times
2N$ matrix $B^{*}$, written as

$B^{*}=\begin{bmatrix}I+\left(A-\Delta\right)+\left(\Delta-I\right)&-\left(\Delta-I\right)\\\
I&O\end{bmatrix}$

can be considered as a state-space matrix of a process on a network, similar
to our LCP process we in (9), with the last $N$ states storing delayed values
of the first $N$ states. The $2N\times 2N$ matrix $B^{*}$ defines the set of
$N$ second-order difference equations, where the governing equation for the
node $i$ position is

$\displaystyle
x_{i}[k+1]=x_{i}[k]+\sum\limits_{j\in\mathcal{N}_{i}}(x_{j}[k]-x_{i}[k])+\left(d_{i}-1\right)\cdot\left(x_{i}[k]-x_{i}[k-1]\right)$

We recognize the second term in (4.4) as an attraction force between
neighbouring nodes with uniform intensity, while in our LCP (9) the attraction
force intensity is proportional to the number of neighbours two adjacent nodes
share. Further, while we propose a repulsive force between adjacent nodes in
(9), node $i$ in (4.4) is repulsed from its previous position $x_{i}[k]$ in
direction of the last position change $(x_{i}[k]-x_{i}[k-1])$.

We implement the weighted intensity of the attractive force as in (9),
ignoring the repulsive force by letting $\delta=0$, and define the $2N\times
2N$ matrix $W^{*}$, corresponding to $B^{*}$,

$W^{*}=\begin{bmatrix}I+\alpha\cdot\Big{(}A\circ
A^{2}+A-\text{diag}\left(\left(A\circ A^{2}+A\right)\cdot
u\right)\Big{)}+\left(\Delta-I\right)&-\left(\Delta-I\right)\\\
I&O\end{bmatrix}.$ (28)

We estimate the number of clusters $c$ in a network from $W^{*}$ similarly as
in the non-back tracking method in Sec. A.3 by counting the number of
eigenvalues in $W^{*}$ with real component larger than
$\sqrt{\lambda_{1}(W^{*})}$. This approach is denoted as LCPc in Section 6.

##  5 Reducing intensity of forces between clusters

The idea behind a group of methods in community detection, called divisive
algorithms, consists of determining the links between nodes from different
clusters. Once these links have been identified, they are removed and thus
only the intra-community links remain [26]. We invoke a similar idea to our
linear clustering process.

An outstanding property of our approach is that the LCP defines the nodal
position as a metric, allowing us to perform clustering in multiple ways. The
position distance between any two, not necessarily adjacent nodes indicates
how likely the two nodes belong to the same cluster. Then, the position metric
also allows us to classify links as either intra- or inter-community. Thus, we
iterate the linear clustering process (9) and, in each iteration, we identify
and scale the weights of the inter-community links.

The attraction and repulsive forces are defined as linear functions of the
position difference between two neighbouring nodes, as presented in Figure 1.
While linear functions greatly simplify the complexity and enable a rigorous
analysis, the linearity of forces introduces some difficulties in the process.
Firstly, as two adjacent nodes are further away, both the attractive and the
repulsive force between them increase in intensity. Similarly, as the
neighbouring nodes are closer on a line, both forces decrease in intensity and
converge to zero as the nodes converge to the same position. Secondly, the
attractive force between any two neighbouring nodes is always of higher
intensity than the repulsive force, causing the process to converge towards
the trivial steady-state.

Non-linearity in the forces can be introduced in the proposed linear
clustering process iteratively by scaling the weights of inter-community links
between iterations, that artificially decreases the strength of forces between
the two nodes from different clusters. In other words, we reduce the
importance of links between nodes from different clusters, based on the
partition from previous iteration.

###  5.1 Scaling the weights of inter-community links

The difference $|\left(y_{2}\right)_{i}-\left(y_{2}\right)_{j}|$ in the
eigenvector $y_{2}$ components of nodes $i$ and $j$ indicates how similar
neighbourhoods of these nodes are. A normalized measure for the difference in
neighbouring nodes $i$ and $j$ is the difference $\left(|r_{i}-r_{j}|\right)$
of their rankings in the sorted eigenvector $\hat{y_{2}}$. Thus, links that
connect nodes with the highest ranking difference are most likely inter-
community links. We define the $N\times N$ scaling matrix $S$ as follows:

$s_{ij}=\begin{cases}1,&\text{if }\lvert r_{j}-r_{i}\rvert<\theta_{r}\\\
\gamma,&\text{otherwise },\end{cases}$ (29)

where the $ij$-th element equals 1 if the absolute value of the ranking
difference between nodes $i$ and $j$ is below a threshold $\theta_{r}$,
otherwise some positive value $0\leq\gamma\leq 1$. Based on the $N\times N$
scaling matrix $S$ in (29), we update the governing equation as follows:

$x[k+1]=\left(I+\tilde{W}-\text{diag}\left(\tilde{W}\cdot u\right)\right)\cdot
x[k],$

where $\tilde{W}=S\circ W$.

Figure 6: Adjacency matrix $A$ of an SSBM network of $N=1000$ nodes, $c=5$
clusters of equal size, with parameters $b_{in}=26$ and $b_{out}=2.25$ (top-
left). Following 4 subfigures present the relabeled adjacency matrix based on
the ranking vector $r$ in iterations $1,5,7$ and $12$, respectively. In each
iteration, the weights of $2\%$ links are scaled (red colour). The weight of
each link is allowed to be scaled once. The relabelled adjacency matrix
$\hat{A}$ after $15$ iterations of scaling weights of links between clusters
(bottom-right).

Scaling the link weights in (29) only impacts the clustering process in (11),
as defined in the equation above. However, modularity-based community
detection, explained in Section 4.2, operates on the $N\times N$ adjacency
matrix $A$ in each iteration. Therefore, our implementation of scaling the
weights of inter-community connections in network helps the process to better
distinguish between clusters (i.e. eventually provides better relabeling in
(20)), without modifying the $N\times N$ adjacency matrix $A$ and, hence,
without negatively affecting the modularity $m$ optimisation in Algorithm 1.
An example of removing links (i.e. $\gamma=0$) is depicted on Figure 6, where
in each iteration weights of $\frac{15}{4}\%$ identified inter-cluster links
are scaled. Scaling the weights of links between clusters significantly
improves the quality of the identified graph partition.

##  6 Benchmarking LCP with other clustering methods

In this section, we benchmark the linear clustering process (9) against
popular clustering algorithms in Appendix A on random networks.

###  6.1 Complexity of LCP

The computational complexity of LCP consists of three parts: the computation
of (i) the $N\times N$ matrix $W$ in (12), (ii) the $N\times 1$ eigenvector
$y_{2}$ of the matrix $W-$diag$(Wu)$ and (iii) the identification of the
clusters based on the sorted eigenvector $\hat{y_{2}}$.

####  6.1.1 Computing the $N\times N$ matrix $W$

The $N\times N$ matrix $A\circ A^{2}$ in (12) requires the highest
computational effort. Generally, computing the square of a matrix involves
$O(N^{3})$ elementary operation, but the zero-one structure of the adjacency
matrix significantly reduces the operations. We provide below an efficient
algorithm for the computation of $A\circ A^{2}$, whose entries determine the
number of $2$-hop walks between any two direct neighbours in the network.

Algorithm 2 Computation of the $N\times N$ matrix $A\circ A^{2}$ .

1:$A$ denotes the adjacency matrix, $N$ denotes number of links, while the set
of node $i$ neighbours is denoted by $\mathcal{N}_{i}$.

2:$A_{s}\leftarrow O_{N\times N}$

3:for $i\leftarrow 1\textrm{ to }N$ do

4: for $j\leftarrow\mathcal{N}_{i}$ do

5: for
$m\leftarrow\left(\mathcal{N}_{j}\setminus\\{1,2,\dots,i\\}\right)\cap\mathcal{N}_{i}$
do

6: $\left(A_{s}\right)_{i,m}\leftarrow\left(A_{s}\right)_{i,m}+1$
$\triangleright$ Account for the 2-hop walk $(i\rightarrow j\rightarrow m)$

7: end for

8: end for

9:end for

10:$A_{s}\leftarrow A_{s}+A_{s}^{T}$

11:return $A_{s}$

We initialize the $N\times N$ matrix $A\circ A^{2}$ with zeros and only
compute elements above the main diagonal, because $A\circ A^{2}$ is symmetric.
The algorithm identifies all 2-hop walks between any two direct neighbours and
accordingly updates the matrix. Let us consider a node $i$ with $d_{i}$
neighbours, denoted as $\mathcal{N}_{i}$. For a neighbouring node
$j\in\mathcal{N}_{i}$, we increment the elements $\left(A\circ
A^{2}\right)_{im}$ by 1, where
$m\in\left(\mathcal{N}_{j}\setminus\\{1,2,\dots,i\\}\right)\cap\mathcal{N}_{i}$,
accounting for 2-hop walks $i\rightarrow j\rightarrow m$. By repeating the
procedure for each node, we compute all the elements above the main diagonal.
Finally, we sum the generated matrix with its transpose to obtain $A\circ
A^{2}$. Since the algorithm 2 is based on incrementing the matrix entries per
each 2-hop walk between direct neighbours, the number of operations equals the
sum $s=u^{T}\cdot\left(A\circ A^{2}\right)\cdot u$ of all elements of $A\circ
A^{2}$

$s=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{N}\lambda_{i}\cdot\lambda_{j}^{2}\cdot
u^{T}\left(x_{i}\circ x_{j}\right)\cdot\left(x_{i}\circ x_{j}\right)^{T}u$
(30)

The eigenvectors of the adjacency matrix $A$ are orthogonal. Therefore
$\left(x_{i}\circ x_{j}\right)^{T}\cdot u=x_{i}^{T}\cdot x_{j}=0$ if $i\neq
j$, otherwise it equals 1 and (30) further simplifies to

$s=\sum\limits_{i=1}^{N}\lambda_{i}^{3},$ (31)

which equals 6 times number of triangles in the network [18, p. 31], because a
2-hop walk between adjacent nodes $i$ and $j$ over a common neighbour $m$ is
equivalent to a triangle $i\rightarrow m\rightarrow j\rightarrow i$. The
computational complexity of $A\circ A^{2}$ thus reduces to $O\left(d_{av}\cdot
L\right)$, as presented in Figure 7. For a given matrix $A\circ A^{2}$, the
computational complexity of the $N\times N$ matrix $W$ is $O(L)$, because (12)
can be transformed into Hadamard product terms (i.e. element-based
operations).

Figure 7: Sum of the cubed eigenvalues $\lambda$ of the adjacency matrix $A$
(blue circles) and product of the average degree $d_{av}$ and number of links
$L$ (red line), for an Erdős–Rényi random graph with $N=300$ nodes, versus the
relative mean degree $\frac{d_{av}}{N}$ (left-figure) and $d_{av}\cdot L$
(right-figure).

####  6.1.2 Computing the $N\times 1$ eigenvector $y_{2}$

The eigenvector $y_{2}$ corresponds to the second largest eigenvalue
$\beta_{2}$ of the $N\times N$ matrix $W-\text{diag}(W\cdot u)$. The largest
eigenvalue $\beta_{1}=1$ corresponds to the eigenvector
$y_{1}=\frac{1}{\sqrt{N}}u$. Computing the eigenvector $y_{2}$ is equivalent
to computing the largest eigenvector of the matrix $W-\text{diag}(W\cdot
u)-\frac{1}{N}\cdot u\cdot u^{T}$, which can be executed using the power
method [18], for a given matrix $W$, with computational complexity
$O\left(L\right)$.

####  6.1.3 Computing the cluster membership function

We apply the recursive algorithm 1 to identify communities based on the
$N\times 1$ eigenvector $y_{2}$. The number of iterations of the algorithm
ideally equals $T=\log_{2}c$, while in worst case scenario there are $c$
iterations. Given a fixed number $c$ of communities, the computational
complexity within an iteration is $O\left(L\right)$, as shown in pseudocode 1.
The number of clusters $c$ may depend upon $N$ and is in worst case equal to
$N$. Thus, computational complexity increases in worst case to $O(N\cdot L)$.

####  6.1.4 Scaling the inter-community links

Between two iterations of the linear clustering process, we identify inter-
community links and scale their weights, as defined in (29). The computational
complexity of this step is $O\left(L\right)$, as the ranking difference of
neighbouring nodes is computer over each link.

Finally, computational complexity of the entire proposed clustering process
equals $O(N\cdot L)$, because $d_{av}=O(N)$.

###  6.2 Clustering performances on stochastic block generated graphs

Figure 8:  The estimated number of clusters (left-figure) in SSBM graphs with
$N=1000$ nodes, average degree $d_{av}=7$ and $c=2$ clusters, for different
values of parameters $b_{in}$ and $b_{out}$. The modularity of the estimated
partitions is presented in right figure. The vertical dashed line indicates
the clustering detectability threshold.

We compare the clustering performance of our LCP with that of clustering
methods introduced in Appendix A, on a same graph generated by the symmetric
stochastic block model (SSBM) with clusters of equal size. All graphs have
$N=1000$ nodes. We vary the parameters $b_{in}$ and $b_{out}$ using (7) in a
way to keep the average degree $d_{av}=7$ fixed. For each SSBM network, we
execute the clustering methods $10^{2}$ times and present the mean number of
estimated clusters and mean modularity of produced partitions in Figures
(8-10). While the non-back tracking algorithm (Sec A.3) and our LCPc (Sec 4.4)
estimate only number of clusters, Newman’s method (Sec A.2), the Louvain
method (Sec A.1) and our LCP (Sec 4.2) and LCPn (Sec 4.3) estimate both number
of clusters and the cluster membership of each node. The attractive strength
$\alpha=0.95$ and the repulsive strength $\delta=10^{-3}$ are used in all
simulations. Weights of $60\%$ links in total are scaled using (29), evenly
over $30$ iterations, where in $i$-th iteration scaled weight is
$\frac{0.05\cdot i}{30}$.

The clustering performance on an SSBM graph with $c=2$ clusters is presented
in Figure 8. The non-back tracking algorithm and our LCPh achieve the best
performance in estimating the number of communities $c$, as shown in left-part
of Figure 8. Further, our LCP outperforms both modularity-based methods, in
identifying the number of communities as well as in modularity $m$.

Figure 8 illustrates a significant difference in performance between our LCP
and the non-back tracking matrix method. The modularity $m$ of our LCP is
plotted in blue, while LCPc, when the number $c$ of communities is known, is
presented in magenta. Our LCP and the other two modularity-based methods
perform poorly in recognising number $c$ of clusters for a wide range of
$b_{in}-b_{out}$, because modularity-based methods generate partitions of
higher modularity than that of the original network (in black colour), but
with different number of communities!

Our LCPn shows that even for a given number of clusters $c$, there exists a
graph partition with a higher modularity than the modularity of the original
network, within theoretically detectable regime.

Figure 9:  The estimated number of clusters (left-figure) in SSBM graphs with
$N=999$ nodes, average degree $d_{av}=7$ and $c=3$ clusters, for different
values of parameters $b_{in}$ and $b_{out}$. The modularity of the estimated
partitions is presented in right figure. The vertical dashed line indicates
the clustering detectability threshold.

In a graph with $c=3$ communities, Figure 9 shows that the non-back tracking
matrix method again outperforms other three methods in identifying number of
communities. Our LCP produces better partitions over the entire range
$b_{in}-b_{out}$ than Newman’s and the Louvain method. As in the case with
$c=2$ clusters, by optimising modularity, the exact number of communities in
SSBM graph cannot always be recovered, because there exist other partitions
with higher number $c$ of clusters with a higher modularity.

Figure 10:  The estimated number of clusters (left-figure) in SSBM graphs with
$N=1000$ nodes, average degree $d_{av}=7$ and $c=8$ clusters, for different
values of parameters $b_{in}$ and $b_{out}$. The modularity of the estimated
partitions is presented in right figure. The vertical dashed line indicates
the clustering detectability threshold.

Finally, Figure 10 shows a similar view for an SSBM graph with $c=8$ clusters.
Our LCP again outperforms the other three methods in estimating number of
communities over the entire range of $b_{in}-b_{out}$ values. In addition, the
modularity of our LCP is superior to the other two methods.

##  7 Conclusion

In this paper, we propose a linear clustering process (LCP) on a network
consisting of an attraction and repulsion process between neighbouring nodes,
proportional to how similar or different their neighbours are. Based on nodal
positions, we are able to estimate both the number $c$ and the nodal
membership of communities. Our LCP outperforms modularity-based clustering
algorithms, such as Newman’s and the Louvain method, while being of the same
computational complexity. The proposed LCP allows estimating the number $c$ of
clusters as accurately as the non-back tracking matrix.

The linear clustering process LCP is described by a matrix
$I+W-\text{diag}(W\cdot u)$, which can be regarded as an operator acting on
the position of nodes, comparable to quantum mechanics (QM). In QM, an
operator describes a dynamical action on a set of particles. Since quantum
mechanical operators are linear, the dynamics are exactly computed via
spectral decomposition. In a same vein, our operator $I+W-\text{diag}(W\cdot
u)$ is linear and describes via attraction and repulsion a most likely
ordering of the position of nodes that naturally leads to clusters, via
spectral decomposition, in particular, via the eigenvector $y_{2}$ in Section
3.3.

## Acknowledgements

The authors are grateful to S. Fortunato for useful comments. This research is
part of NExTWORKx, a collaboration between TU Delft and KPN on future
communication networks.

##  Appendix A Clustering algorithms

###  A.1 Louvain Method

The Louvain method is a simple, yet powerful heuristic clustering algorithm,
proposed by Blondel et al. [13]. The method is based on an iterative,
unsupervised, two-step procedure that optimizes modularity $m$. Initially, a
directed graph $G$ with an $N\times N$ weighted adjacency matrix $M$ is
partitioned in $N$ clusters, where each node constitutes its own cluster or
community.

In the first stage, the algorithm examines how the graph modularity $m$
changes if a node $i$ would be assigned to a community of its neighbouring
node $j\in\mathcal{N}_{i}$. The modularity gain $\Delta m$ in case node $i$ is
assigned to community $h$ of adjacent node $j$ has been determined in [13] as

$\displaystyle\Delta
m=\left(\frac{\sum_{\text{in}}+2\sum_{l:C_{lj}=1}M_{il}}{2L}-\left(\frac{\sum_{\text{tot}}+d_{i}}{2L}\right)^{2}\right)-\left(\frac{\sum_{\text{in}}}{2L}-\left(\frac{\sum_{\text{tot}}}{2L}\right)^{2}-\left(\frac{d_{i}}{2L}\right)^{2}\right),$
(32)

where the sum of the weights of intra-community links in $h$ is
$\sum_{\text{in}}$, while $\sum_{\text{tot}}$ denotes the sum of the weights
of all links in $G$ incident to any node in community $h$. Node $i$ is
assigned to the community with the largest positive gain in modularity $m$. In
case all computed gains $\Delta m$ are either negative or smaller than a
predefined small positive threshold value, node $i$ remains in its original
community. The first stage ends when modularity $m$ cannot be further
increased by re-assigning nodes to communities of neighbours.

In the second stage of an iteration, the weighted graph from the first stage
is transformed into a new weighted graph, where each community is presented by
a node. The link weight between two nodes $h$ and $g$ equals the sum of
weights of all links between communities $h$ and $g$ in the graph from the
first stage. Furthermore, the weight of a self-loop of node $g$ in the new
graph equals the sum of weights of all intra-community links in cluster $g$ of
the graph from the previous stage. The new graph is provided to the first
stage in the next iteration. The algorithm stops when modularity $m$ cannot be
increased further. The time complexity of the Louvain method is linear in the
number of links $O(L)$ on sparse graphs [13].

###  A.2 Newman’s Method of Optimal Modularity

Newman [12] proposed a clustering algorithm that is based on modularity
optimisation. The algorithm starts with estimating the bisection of a graph
$G$, generating the highest modularity $m$ from (3), that can be rewritten as
follows:

$m=\frac{1}{4L}y^{T}\cdot M\cdot y,$ (33)

where the $N\times 1$ vector $y$ is composed of values $1$ and $-1$, denoting
cluster membership of each node, while the $N\times N$ modularity matrix
$M=A-\frac{1}{2L}\cdot d\cdot d^{T}$ has the following eigenvalue
decomposition

$M=\sum\limits_{i=1}^{N}\zeta_{i}\cdot z_{i}\cdot z_{i}^{T},$ (34)

where the $N\times 1$ eigenvector $z_{i}$ corresponds to the $i$-th eigenvalue
$\zeta_{i}$. Further, the vector $y=\sum_{j=1}^{N}(z_{j}^{T}\cdot y)\cdot
z_{j}$ can be written as a linear combination of eigenvectors
$\\{z_{i}\\}_{1\leq i\leq N}$, which transforms (33) to

$m=\frac{1}{4L}\sum\limits_{i=1}^{N}\zeta_{i}\cdot(z_{j}^{T}\cdot y)^{2}.$
(35)

In order to maximise the modularity $m$, Newman [12] proposed to define
$y_{i}=1$, in case $(z_{1})_{i}>0$, otherwise $y_{i}=-1$. In a next iteration,
the same procedure of spectral division into two partitions is repeated on
both sub-graphs. However, using only the block sub-matrix of $M$,
corresponding to cluster $g$ in next iteration would not take into account
inter-community links. Instead, for the estimated cluster $g$, the modularity
matrix $M_{g}$ is updated as

$M_{g}=m_{ij}-\left(\sum_{k\in g}m_{ik}\right)\cdot\delta_{ij},$ (36)

where Kronecker delta $\delta_{ij}=1$ if $i=j$, otherwise $\delta_{ij}=0$. The
algorithm stops when the modularity $m$ cannot be further improved.

###  A.3 Non-back tracking matrix

The non-back tracking clustering method estimates the number of clusters in a
network, based on the spectrum of the non-back tracking matrix $B$, that
contains information about $2$-hop directed walks in a network $G$, that are
not closed [17]. Given an undirected network $G(\mathcal{N},\mathcal{L})$, for
each link $i\sim j$ between nodes $i$ and $j$, two directed links
$(i\rightarrow j)$ and $(j\rightarrow i)$ are created. By transforming each
link in $G$ into a bi-directional link pair, we compose in total $2L$ links.
The $2L\times 2L$ non-back tracking matrix $B$ is defined as follows:

$B_{(u\rightarrow v),(w\rightarrow z)}=\begin{cases}1&\text{if $v=w$ and
$u\neq z$}\\\ 0&\text{otherwise},\end{cases}$ (37)

where $v,w,z\in\mathcal{N}$. Since the non-back tracking matrix $B$ is
asymmetric, its eigenvalues are generally complex. Furthermore, a vast
majority of eigenvalues lie within a circle in complex plain, with centre at
the origin and with radius equal to the square root of the largest eigenvalue.
Krzakala et al. [17] hypothesized that the number of clusters in $G$ equals
the number of real-valued eigenvalues outside the circle. Computing the
eigenvalues of the non-back tracking matrix $B$ is of computational complexity
$O(L^{3})$. However, the complexity can be reduced to $O(N^{3})$, as explained
in [10, p. 20]. The non-back tracking matrix method is denoted as NBTM in
Section 6.

##  Appendix B LCP in continuous time

We explain the physical intuition of our clustering process in continuous time
$t$, where the position $x_{i}(t)$ of a node $i$ is assumed to change
continuously with time $t$. The change $x_{i}(t+\Delta t)-x_{i}(t)$ in
position of node $i$ at time $t$ for small increments $\Delta t$ is
proportional to the sum over neighbours $j$ of the difference
$x_{j}(t)-x_{i}(t)$ in position weighted by the resultant force between
attraction and repulsion:

$\frac{dx_{i}(t)}{dt}=\sum\limits_{j\in\mathcal{N}_{i}}\bigg{(}\frac{\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}{d_{j}d_{i}}-\frac{\frac{1}{2}\cdot\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{d_{j}d_{i}}\bigg{)}\cdot\Big{(}x_{j}(t)-x_{i}(t)\Big{)}$
(38)

where $\alpha$ and $\delta$ are, in the continuous-time setting, the rates
(with units s-1) for attraction and repulsion, respectively. The law (38) of
the nodal positioning $x_{i}\left(t\right)$ for each $i\in\mathcal{N}$
deviates from physical repulsion between charged particles, where the force is
proportional to $\left(x_{j}(t)-x_{i}(t)\right)^{-b}$ for some positive number
$b$. The important advantage of the law (38) is its linearity that allows an
exact mathematical treatment. The linear dynamic process (38) is proportional
to the underlying graph, which we aim to cluster; a non-linear law depends
intricately on the underlying graph and may result in a lesser clustering. The
drawback of the linear dynamic process (38), as investigated below in Section
3.3, lies in the steady state, where the attractive and repulsive forces are
precisely in balance.

After dividing both sides by $\delta$,

$\frac{dx_{i}(t)}{d\left(\delta
t\right)}=\sum\limits_{j\in\mathcal{N}_{i}}\bigg{(}\frac{\frac{\alpha}{\delta}\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}{d_{i}d_{j}}-\frac{\frac{1}{2}\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{d_{i}d_{j}}\bigg{)}\cdot\Big{(}x_{j}(t)-x_{i}(t)\Big{)}$

and defining the normalized time by $t^{\ast}=\delta t$ and the effective
attraction rate $\tau=\frac{\alpha}{\delta}$, the governing equation (38)
reduces to

$\frac{dx_{i}(t^{\ast})}{dt^{\ast}}=\sum\limits_{j\in\mathcal{N}_{i}}\Bigg{(}\frac{\tau\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}{d_{i}d_{j}}-\frac{\left(\frac{\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|}{2}-1\right)}{d_{i}d_{j}}\Bigg{)}\cdot\Big{(}x_{j}(t^{\ast})-x_{i}(t^{\ast})\Big{)}$
(39)

The position $x_{i}(t^{\ast})$ of node $i$ is now expressed in the
dimensionless time $t^{\ast}$, where the actual time
$t=\frac{t^{\ast}}{\delta}$ is measured in units of $\frac{1}{\delta}$. By
scaling or normalizing the time, the repulsion strength or rate $\delta$ has
been eliminated, illustrating that the clustering process only depends upon
one parameter, the effective attraction rate $\tau$. Relation (1) indicates
that the weight of the position difference

$w_{ij}=\frac{\tau\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)-\left(\frac{\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|}{2}-1\right)}{d_{i}d_{j}}$

lies in the interval
$\left(-\frac{\frac{d_{i}+d_{j}}{2}-1}{d_{i}d_{j}},\frac{\tau}{d_{i}}\right)$
and that the elements $w_{ij}=w_{ji}$ define the symmetric $N\times N$ weight
matrix $W$, which is specified in (12) below. Although symmetry is physically
not required 111The process described by
$\frac{dx_{i}(t)}{dt}=\sum\limits_{j\in\mathcal{N}_{i}}\frac{\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)-\delta\cdot\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|-1\right)}{d_{j}d_{i}}\cdot\Big{(}x_{j}(t)-x_{i}(t)\Big{)}$
also works., the analysis below is greatly simplified, because eigenvalues and
eigenvectors of a symmetric matrix are real.

We rewrite the law (39) as

$\displaystyle\frac{dx_{i}(t^{\ast})}{dt^{\ast}}=\sum\limits_{j\in\mathcal{N}_{i}}w_{ij}x_{j}(t^{\ast})-x_{i}(t^{\ast})\sum\limits_{j\in\mathcal{N}_{i}}w_{ij}=\sum_{j=1}^{N}a_{ij}w_{ij}x_{j}(t^{\ast})-x_{i}(t^{\ast})v_{i}$

where
$v_{i}=\sum\limits_{j\in\mathcal{N}_{i}}w_{j}=\sum_{j=1}^{N}a_{ij}w_{ij}$ is
independent of time $t^{\ast}$. The cluster positioning law (38) for the
vector $x(t^{\ast})$ in continuous time is, in matrix form,

$\frac{dx(t^{\ast})}{dt^{\ast}}=(A\circ W-\text{diag}(v))x(t^{\ast})$ (40)

where the Hadamard product [21] is denoted by $\circ$ and the vector
$v=(A\circ W)u$. The corresponding solution of (40) is [27, eq. (6)]

$x(t^{\ast})=e^{(A\circ W-\text{diag}((A\circ W)u))t^{\ast}}x(0)$ (41)

which illustrates that a steady state is reached, provided that the real part
of the largest eigenvalue of the matrix $H=(A\circ W-$diag$((A\circ W)u))$ is
not positive.

##  Appendix C Proof of Theorems

###  C.1 Proof of Theorem 1

Similarly as in Section B, we rewrite the sum over all neighbours in the
governing equation (9) in terms of the elements of the $N\times N$ adjacency
matrix $A$:

$x_{i}[k+1]-x_{i}[k]=\sum\limits_{j=1}^{N}\frac{a_{ij}}{d_{i}d_{j}}\Big{(}x_{j}[k]-x_{i}[k]\Big{)}\Big{(}\alpha\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|-\frac{1}{2}\delta\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|\right)\Big{)}.$
(42)

Firstly, we denote the $N\times 1$ vector $\tilde{d}=\Delta^{-1}\cdot u$
composed of the inverse nodal degrees:

$\tilde{d}=\begin{bmatrix}\frac{1}{d_{1}}&\frac{1}{d_{2}}&\dots&\frac{1}{d_{N}}\end{bmatrix}^{T}$
(43)

In the sequel, we will deduce the corresponding matrix form of (42). With (1)
and (2), the degree $d_{i}$ of node $i$ distracted by the number of common
neighbours between nodes $i$ and $j$ (2), equals the number of node $i$
neighbours, not adjacent to node $j$:

$\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|=\left(d\cdot
u^{T}-A^{2}\right)_{ij}.$ (44)

Similarly, the number of node $j$ neighbours that do not share link with node
$i$ has following matrix form:

$\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|=\left(u\cdot
d^{T}-A^{2}\right)_{ij}.$ (45)

Finally, the position difference $(x_{j}[k]-x_{i}[k])$ between nodes $i$ and
$j$ at time $k$ equals the $ij$-th element of the matrix below:

$\left(x_{j}[k]-x_{i}[k]\right)=\left(u\cdot x^{T}[k]-x[k]\cdot
u^{T}\right)_{ij},$ (46)

while dividing by node $i$ ($j$) degree $d_{i}$ ($d_{j}$) is equivalent to
product with the $ij$-th element of the $N\times N$ matrix $(\tilde{d}\cdot
u^{T})_{ij}$ and $(u\cdot\tilde{d}^{T})_{ij}$, respectively. By implementing
matrix notations (2), (44), (45) and (46) into the governing equation (42) and
by applying the distributive property of the Hadamard product [21, p. 477] we
obtain:

$\displaystyle x[k+1]-x[k]=$ $\displaystyle\bigg{(}\left(u\cdot
x^{T}[k]-x[k]\cdot u^{T}\right)\circ A\circ$ (47)
$\displaystyle\left(u\cdot\tilde{d}^{T}\right)\circ\left(\tilde{d}\cdot
u^{T}\right)\circ\Big{(}\left(\alpha+\delta\right)\cdot\left(A^{2}+A\right)-$
$\displaystyle\frac{1}{2}\delta\cdot\left(u\cdot d^{T}+d\cdot
u^{T}\right)\Big{)}\bigg{)}\cdot u.$

We define the $N\times N$ topology-based matrix $W$ as follows:

$W=A\circ\left(u\cdot\tilde{d}^{T}\right)\circ\left(\tilde{d}\cdot
u^{T}\right)\circ\Big{(}\left(\alpha+\delta\right)\cdot\left(A^{2}+A\right)-\frac{1}{2}\delta\left(u\cdot
d^{T}+d\cdot u^{T}\right)\Big{)}.$ (48)

Using the distributive property of a Hadamard product [21, p. 477], we develop
the equation (48) further:

$\begin{split}W=&\left(\alpha+\delta\right)\cdot\Big{(}\left(u\cdot\tilde{d}^{T}\right)\circ\left(\tilde{d}\cdot
u^{T}\right)\circ A\circ\left(A^{2}+A\right)\Big{)}-\\\
&\frac{1}{2}\delta\Big{(}A\circ\left(u\cdot\tilde{d}^{T}\right)\circ\left(\tilde{d}\cdot
u^{T}\right)\circ\left(u\cdot d^{T}\right)\Big{)}-\\\
&\frac{1}{2}\delta\Big{(}A\circ\left(u\cdot\tilde{d}^{T}\right)\circ\left(\tilde{d}\cdot
u^{T}\right)\circ\left(d\cdot u^{T}\right)\Big{)}.\end{split}$ (49)

Since the Hadamard product is commutative [21, p. 477], we can reorder the
products in previous equation. The Hadamard product
$\left(u\cdot\tilde{d}^{T}\right)\circ\left(u\cdot d^{T}\right)$ equals all-
one matrix $J$. Similarly, the product $\left(\tilde{d}\cdot
u^{T}\right)\circ\left(d\cdot u^{T}\right)=J$ . We further transform the
Hadamard product of $\left(A\circ A^{2}+A\right)$ and the outer products
$\left(u\cdot\tilde{d}^{T}\right)$ and $\left(\tilde{d}\cdot u^{T}\right)$
into product with diagonal matrices $\Delta^{-1}\cdot\left(A\circ
A^{2}+A\right)\cdot\Delta^{-1}$. Thus, equation (49) transforms to (12).
Substituting (48) into (47) yields

$x[k+1]-x[k]=\bigg{(}\left(u\cdot x^{T}[k]-x[k]\cdot u^{T}\right)\circ
W\bigg{)}\cdot u.$ (50)

The Hadamard product of a matrix with an outer product of two vectors is
equivalent to the product with diagonal matrices of vectors composing the
outer product [21, p. 477]. Thus, we further transform the governing equation
(50):

$x[k+1]-x[k]=W\cdot\text{diag}\left(x[k]\right)\cdot
u-\text{diag}\left(x[k]\right)\cdot\left(W\cdot u\right),$ (51)

where the last term $\text{diag}\left(x[k]\right)\cdot\left(W\cdot u\right)$
represents the Hadamard product of two vectors, $x[k]\circ\left(W\cdot
u\right)$ and can be presented as $\text{diag}\left(W\cdot u\right)\cdot
x[k]$. Thus, the equation transforms into (11) which completes the proof.
$\hfill\square$

###  C.2  Proof of Property 2

We observe that

$W\cdot u-\text{diag}\left(W\cdot u\right)\cdot u=0$

implying that the all-one vector $u$ is an eigenvector of the matrix
$W-\text{diag}\left(W\cdot u\right)$ belonging to the zero eigenvalue.
Therefore, the $N\times N$ matrix $I+W-\text{diag}\left(W\cdot u\right)$ has
an eigenvalue $1$ corresponding to the all-one vector $u$.

By the Perron-Frobenius theorem [28] for a non-negative matrix, the principal
eigenvector, belonging to the largest eigenvalue, has non-negative components.
Since the eigenvector $u$ has non-negative components and all eigenvectors of
a symmetric matrix are orthogonal, it follows that the all-one vector $u$ is
the Perron or principal eigenvector belonging to the largest eigenvalue 1 of
the matrix $I+W-\text{diag}\left(W\cdot u\right)$ and, thus, all other real
eigenvalues are, in absolute value, smaller than 1. $\hfill\square$

###  C.3 Proof of Property 3

The non-negativity of the matrix $I+W-\text{diag}\left(W\cdot u\right)$
implies that $w_{ij}\geq 0$ for $i\neq j$ and
$1+w_{ii}-\sum_{k=1}^{N}w_{ik}\geq 0$, hence,

$1\geq\sum_{k=1;k\neq i}^{N}w_{ik}\geq 0$

Equivalently, the symmetric matrix $W-\text{diag}\left(W\cdot u\right)$ has
positive off-diagonal elements, but negative diagonal elements, similarly to
the infinitesimal generator of a Markov chain (which is minus a weighted
Laplacian [29]). Introducing the explicit expression (13) and requiring that
each element $w_{ij}$ is non-negative,

$w_{ij}=a_{ij}\frac{\alpha\cdot\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)-\delta\cdot\left(\frac{\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|}{2}-1\right)}{d_{i}d_{j}}\geq
0$

leads to

$\frac{\alpha}{\delta}\geq\frac{1}{2}\cdot\frac{\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}$

which holds for any $i,j\neq i\in\mathcal{N}$. With (1), (2) and
$d_{i}=(Au)_{i}=(A^{2})_{ii}$, the condition for the ratio
$\frac{\alpha}{\delta}$ becomes

$\frac{\alpha}{\delta}\geq\max_{i,j\neq
i\in\mathcal{N}}\frac{1}{2}\cdot\frac{\left(\left|\mathcal{N}_{j}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{j}\right|-2\right)}{\left(\left|\mathcal{N}_{j}\cap\mathcal{N}_{i}\right|+1\right)}=\max_{i,j\neq
i\in\mathcal{N}}\frac{d_{i}+d_{j}}{2\left((A^{2})_{ij}+1\right)}-1$

which simplifies to

$\frac{\alpha}{\delta}\geq d_{\max}-1$ (52)

We write $\sum_{k=1;k\neq i}^{N}w_{ik}$ with (13) as

$\sum_{k=1;k\neq
i}^{N}w_{ik}=\frac{1}{d_{i}}\sum_{k=1}^{N}a_{ik}\bigg{(}\frac{\alpha\cdot\left(\left|\mathcal{N}_{k}\cap\mathcal{N}_{i}\right|+1\right)}{d_{k}}-\frac{\frac{\delta}{2}\cdot\left(\left|\mathcal{N}_{k}\setminus\mathcal{N}_{i}\right|+\left|\mathcal{N}_{i}\setminus\mathcal{N}_{k}\right|-2\right)}{d_{k}}\bigg{)}$

Introducing (1) and (2),

$\displaystyle\sum_{k=1;k\neq i}^{N}w_{ik}$
$\displaystyle=\frac{1}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\alpha\cdot\left((A^{2})_{ik}+1\right)-\frac{\delta}{2}\cdot\left(d_{i}+d_{k}-2\left((A^{2})_{ik}+1\right)\right)}{d_{k}}$
$\displaystyle=\frac{\alpha}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}-\frac{\delta}{2d_{i}}\sum_{k=1}^{N}a_{ik}\left(\frac{d_{i}}{d_{k}}+1-\frac{2\left((A^{2})_{ik}+1\right)}{d_{k}}\right)$

leads, with $d_{i}=\sum_{k=1}^{N}a_{ik}$, to

$\displaystyle\sum_{k=1;k\neq
i}^{N}w_{ik}=\frac{\alpha+\delta}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}-\frac{\delta}{2}-\frac{\delta}{2}\sum_{k=1}^{N}\frac{a_{ik}}{d_{k}}$

The second condition $\sum_{k=1;k\neq i}^{N}w_{ik}\leq 1$,

$\frac{\alpha+\delta}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}-\frac{\delta}{2}-\frac{\delta}{2}\sum_{k=1}^{N}\frac{a_{ik}}{d_{k}}\leq
1$

must hold for all $i\in\mathcal{N}$, which translates to

$\displaystyle 1$
$\displaystyle\geq\max_{i\in\mathcal{N}}\left(\frac{\alpha+\delta}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}-\frac{\delta}{2}-\frac{\delta}{2}\sum_{k=1}^{N}\frac{a_{ik}}{d_{k}}\right)$
$\displaystyle\geq\left(\alpha+\delta\right)\max_{i\in\mathcal{N}}\frac{1}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}-\frac{\delta}{2}-\frac{\delta}{2}\min_{i\in\mathcal{N}}\sum_{k=1}^{N}\frac{a_{ik}}{d_{k}}$

With $\left((A^{2})_{ik}+1\right)\leq d_{k}$ if $a_{ik}=1$, we have
$\frac{1}{d_{i}}\sum_{k=1}^{N}a_{ik}\frac{\left((A^{2})_{ik}+1\right)}{d_{k}}\leq
1$, while
$\frac{d_{i}}{d_{\min}}\geq\sum_{k=1}^{N}\frac{a_{ik}}{d_{k}}\geq\frac{d_{i}}{d_{\max}}$.
Hence, the second condition becomes

$1\geq\alpha+\frac{\delta}{2}\left(1-\frac{d_{\min}}{d_{\max}}\right)$ (53)

illustrating that $\alpha\leq 1$. Combining the two conditions (52) and (53)
into a linear set of inequalities

$\left\\{\begin{array}[c]{c}0\geq-\alpha+\delta\left(d_{\max}-1\right)\\\
1\geq\alpha+\frac{\delta}{2}\left(1-\frac{d_{\min}}{d_{\max}}\right)\end{array}\right.$

or in matrix form, where $\curlyeqsucc$ denotes componentwise inequalities
[30, p. 32, 40]

$\left[\begin{array}[c]{c}0\\\
1\end{array}\right]\curlyeqsucc\left[\begin{array}[c]{cc}-1&\left(d_{\max}-1\right)\\\
1&\frac{1}{2}\left(1-\frac{d_{\min}}{d_{\max}}\right)\end{array}\right]\left[\begin{array}[c]{c}\alpha\\\
\delta\end{array}\right]$

yields, after inversion, the bounds (17) and (18). $\hfill\square$

## References

  * [1] A.-L. Barabási, “Network science,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 371, no. 1987, p. 20120375, 2013.
  * [2] M. Newman, Networks.  Oxford university press, 2018.
  * [3] B. Karrer, M. E. Newman, and L. Zdeborová, “Percolation on sparse networks,” Physical review letters, vol. 113, no. 20, p. 208702, 2014.
  * [4] B. Prasse and P. Van Mieghem, “The viral state dynamics of the discrete-time nimfa epidemic model,” IEEE Transactions on Network Science and Engineering, vol. 7, no. 3, pp. 1667–1674, 2019.
  * [5] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, “Epidemic processes in complex networks,” Reviews of modern physics, vol. 87, no. 3, p. 925, 2015.
  * [6] A. D. Sánchez, J. M. López, and M. A. Rodriguez, “Nonequilibrium phase transitions in directed small-world networks,” Physical review letters, vol. 88, no. 4, p. 048701, 2002.
  * [7] H. E. Stanley, Phase transitions and critical phenomena, vol. 7.  Clarendon Press, Oxford, 1971.
  * [8] I. Jokić and P. Van Mieghem, “Linear processes on complex networks,” Journal of Complex Networks, vol. 8, no. 4, p. cnaa030, 2020.
  * [9] S. Fortunato, “Community detection in graphs,” Physics reports, vol. 486, no. 3-5, pp. 75–174, 2010.
  * [10] G. Budel and P. Van Mieghem, “Detecting the number of clusters in a network,” Journal of Complex Networks, vol. 8, no. 6, p. cnaa047, 2020.
  * [11] M. E. Newman, “Communities, modules and large-scale structure in networks,” Nature physics, vol. 8, no. 1, pp. 25–31, 2012.
  * [12] M. E. Newman, “Modularity and community structure in networks,” Proceedings of the national academy of sciences, vol. 103, no. 23, pp. 8577–8582, 2006.
  * [13] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,” Journal of statistical mechanics: theory and experiment, vol. 2008, no. 10, p. P10008, 2008.
  * [14] A. Lancichinetti and S. Fortunato, “Community detection algorithms: a comparative analysis,” Physical review E, vol. 80, no. 5, p. 056117, 2009\.
  * [15] T. P. Peixoto, “Hierarchical block structures and high-resolution model selection in large networks,” Physical Review X, vol. 4, no. 1, p. 011047, 2014.
  * [16] L. Danon, A. Diaz-Guilera, J. Duch, and A. Arenas, “Comparing community structure identification,” Journal of statistical mechanics: Theory and experiment, vol. 2005, no. 09, p. P09008, 2005.
  * [17] F. Krzakala, C. Moore, E. Mossel, J. Neeman, A. Sly, L. Zdeborová, and P. Zhang, “Spectral redemption in clustering sparse networks,” Proceedings of the National Academy of Sciences, vol. 110, no. 52, pp. 20935–20940, 2013.
  * [18] P. Van Mieghem, Graph spectra for complex networks.  Cambridge University Press, 2010.
  * [19] P. Van Mieghem, K. Devriendt, and H. Cetinay, “Pseudoinverse of the laplacian and best spreader node in a network,” Physical Review E, vol. 96, no. 3, p. 032311, 2017.
  * [20] M. E. Newman and M. Girvan, “Finding and evaluating community structure in networks,” Physical review E, vol. 69, no. 2, p. 026113, 2004.
  * [21] R. A. Horn and C. R. Johnson, Matrix analysis.  Cambridge university press, 2012.
  * [22] P. W. Holland, K. B. Laskey, and S. Leinhardt, “Stochastic blockmodels: First steps,” Social networks, vol. 5, no. 2, pp. 109–137, 1983.
  * [23] A. Decelle, F. Krzakala, C. Moore, and L. Zdeborová, “Inference and phase transitions in the detection of modules in sparse networks,” Physical Review Letters, vol. 107, no. 6, p. 065701, 2011.
  * [24] A. Decelle, F. Krzakala, C. Moore, and L. Zdeborová, “Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications,” Physical Review E, vol. 84, no. 6, p. 066106, 2011.
  * [25] O. Angel, J. Friedman, and S. Hoory, “The non-backtracking spectrum of the universal cover of a graph,” Transactions of the American Mathematical Society, vol. 367, no. 6, pp. 4287–4318, 2015.
  * [26] M. Girvan and M. E. Newman, “Community structure in social and biological networks,” Proceedings of the national academy of sciences, vol. 99, no. 12, pp. 7821–7826, 2002.
  * [27] P. Van Mieghem, “Approximate formula and bounds for the time-varying susceptible-infected-susceptible prevalence in networks,” Physical Review E, vol. 93, no. 5, p. 052312, 2016.
  * [28] P. Van Mieghem, X. Ge, P. Schumm, S. Trajanovski, and H. Wang, “Spectral graph analysis of modularity and assortativity,” Physical Review E, vol. 82, no. 5, p. 056113, 2010.
  * [29] P. Van Mieghem, Performance analysis of complex networks and systems.  Cambridge University Press, 2014.
  * [30] S. Boyd, S. P. Boyd, and L. Vandenberghe, Convex optimization.  Cambridge university press, 2004.